№ 4
ИЗВЕСТИЯ АКАДЕМИИ НАУК ЭНЕРГЕТИКА
2008
УДК 621.039.586:536.24
© 2008 г. ЧУДАНОВ В.В., АКСЕНОВА А.Е., ПЕРВИЧКО В.А.
РАЗРАБОТКА МОДЕЛЕЙ И МЕТОДОВ 3D ТЕПЛО- И ГИДРОДИНАМИКИ ДЛЯ МОДЕЛИРОВАНИЯ ДВУХФАЗНЫХ ПУЗЫРЬКОВЫХ И КАПЕЛЬНЫХ ТЕЧЕНИЙ С УЧЕТОМ СИЛ ПОВЕРХНОСТНОГО НАТЯЖЕНИЯ
В статье рассматривается вычислительная методика исследования пузырьковых и капельных течений. Предлагаются эффективные вычислительные CFD алгоритмы для 3D расчетов двухфазных течений с явным выделением межфазной границы и учетом сил поверхностного натяжения. Акцент сделан на валидации и верификации вычислительной математической модели с использованием доступных из литературы численных и экспериментальных тестов.
Введение. Несмотря на длительную историю исследования двухфазных и парожид-костных течений удовлетворительного вычислительного подхода для определения движения разных фаз и поведения двухфазных смесей до сих пор нет. Одной из причин такого положения является отсутствие детальной картины взаимодействия пузырьков с несущей средой и между собой. Эти задачи сложны для прямого теоретического исследования из-за большого числа неопределенных параметров и особенно сложны для нестационарных условий взаимодействия.
В этой связи актуальны развитие вычислительных методов и алгоритмов для решения задач тепло- и массопереноса в несжимаемых средах с учетом сил поверхностного натяжения и разработка вычислительных модулей для трехмерного CFD кода [1], отвечающего современному уровню развития вычислительных средств и численных методов.
Один из часто встречающихся видов двухфазных CFD течений - пузырьковое течение. Этот режим реализуется при развитом кипении жидкости на различно ориентированных поверхностях, двухфазном течении в вертикальных и горизонтальных каналах. Разрабатываемые и используемые методы прямого численного моделирования двухфазных течений нуждаются в тестировании на доступных из литературы экспериментальных и аналитических данных, в частности на тестах, описывающих установившееся движение одиночного газового пузырька и группы пузырьков в большом объеме несжимаемой жидкости. Необходимы эксперименты по описанию всплытия пузырька в каналах различной формы (цилиндрический, щелевой каналы) и эффектов взаимовлияния пузырьков при слиянии (коалесценции) и дроблении. В связи с этим актуальна разработка системы тестов (матрицы валидации) для верификации создаваемого программного обеспечения прямого моделирования пузырьковых течений и обобщение результатов в виде карты режимов течений всплывающего пузырька и группы пузырьков.
В работе приведены особенности разработанной в ИБРАЕ РАН методики моделирования пузырьковых течений и результаты ее валидации на основе экспериментов, опубликованных [2-4]. Особое внимание уделено течениям при малых числах Рейнольдса, а именно экспериментальным данным по всплытию одного и группы пузырьков в плоском канале и экспериментальным данным по взаимодействию капель жидкости с малой разностью плотностей с учетом слияния (коалесценции) и разрыва капель. Кроме того, тестирование разрабатываемого программного обеспечения выполнено на таких из-
вестных тестах как капиллярная волна, неустойчивость Рейлея-Тейлора, задача об обрушении плотины и т.д. [5].
СЕБ вычислительная методика. Данный раздел посвящен новым методам и алгоритмам для решения 3Б тепло- и гидродинамики в несжимаемых жидкостях. Среди них алгоритмы для решения несжимаемой динамики жидкости, монотонные многомерные схемы ТУБ-типа для решения уравнения адвекции, эффективные алгоритмы решения эллиптического уравнения для поправки давления. Для моделирования поведения пара внутри пузырька предлагается модель сжимаемых сред при малых числах Маха [6]. Для моделирования поведения жидкости вне пузырька используются нестационарные уравнения Навье-Стокса в естественных переменных для несжимаемой жидкости вместе с уравнением энергии [7].
Для наблюдения за поверхностью раздела двухфазной системы использованы:
- модифицированный метод набора уровней для прямого численного моделирования двухфазной системы типа жидкость-пар [1];
- монотонная многомерная схема адвекции с малой схемной диффузией с использованием подсеточного моделирования;
- эффективные вычислительные алгоритмы для решения уравнения давления в гидродинамике несжимаемой жидкости с переменной плотностью.
Алгоритмы, методы и программное обеспечение были тестированы при широком наборе тестов [8-10].
Вычислительная методика для решения уравнений Навье-Стокса. Для моделирования поведения жидкости вне пузырьков использованы нестационарные уравнения Навье-Стокса в естественных переменных [7] для несжимаемой жидкости вместе с уравнением энергии
[ йрь/йг = Р + V и+ + С8Б;
I
I Шу и = 0;
^^ + div (puh) = div (к grad T) h = [ с (£) 4. (2)
dt J
0
Перечислим основные особенности развитого численного алгоритма [1, 8], включенные в предложенный метод. Численная реализация схем расщепления для уравнений Навье-Стокса выполнена как предиктор-корректор процедура с коррекцией для поправки давления
[р(и - и )/т] + (C(и)
-divvgrad) и + gradp -CSF =0; (3)
.. f 1 , о j 1 n + 1/2 n + 1 n + 1/2 т.е. ,„4
divh[pgradhspj = Tdivhu ; u =u -pgradh8p. (4)
T-. .S
В случае схемы с итерациями система уравнений с учетом итераций по v имеет вид
1/2
, и
Pol
(v __—uY x + ^f^ + (с(и) -div v grad ) и = P g + CSF;
т /V Po J Po Po
S n + 1/2n , i. \ s л \ /»S n + 1/2N ~s (c\
v - и j Apj Л Apj,. f v - v j ,P (5)
■ — + 1 + — div- = -div v grad—;
l div v = 0.
gradl 1 + — + I 1 + — div I- = -div v grad—;
1 PoJ V PoJ V т J Po
Принимая во внимание последнюю строку в (5), получаем следующее уравнение для коррекции скорости
„S + 1 V
\
/
(6)
Для решения задачи конвекции создана регуляризованная нелинейная монотонная операторная схема расщепления [11]. Специальная аппроксимация конвективных слагаемых С(и) использована для получения дискретного конвективного оператора, который является кососимметричным и не дает вклада в кинетическую энергию (т.е. энергетически нейтрален).
Данная схема обеспечивает второй порядок по пространству и первый по времени. Алгоритм устойчив при достаточно большом шаге интегрирования по времени.
Результаты тестирования метода на двумерных и трехмерных тестах приведены в [7].
Быстрый вычислительный модуль для решения уравнения давления. Для решения эллиптических уравнений с переменными коэффициентами, к примеру, уравнения Пуассона для поправки давления Ър вида
обычно применяют метод сопряженного градиента, который занимает 90% времени центрального процессора на один расчетный шаг. Поэтому важен поиск алгоритмов, требующих меньшего количества времени центрального процессора на один расчетный шаг.
Альтернативой методу сопряженных градиентов может служить итерационный метод Ричардсона с Чебышевским набором параметров, использующий в качестве предо-буславливателя БПФ солвер для оператора Лапласа с постоянным множителем. Применение этого метода для решения эллиптических уравнений с переменными коэффициентами позволяет достичь 50-ти кратного ускорения по сравнению с обычно используемым методом сопряженных градиентов [8].
Кривизна и поверхностное натяжение. Термин капиллярного источника, появляющийся в уравнениях движения может быть аппроксимирован в виде [12]
где n обозначает единичный вектор нормали к поверхности раздела; K - положительная константа; к - локальная кривизна, определяемая как
к = -Vn = -V(Vx/| Vxl).
Интегрирование уравнения по интерфейсной области (ds) дает модель непрерывной поверхностной силы (CSF), см. [13]
|ок5( х - Xf) nds,
где о - поверхностное натяжение, которое предполагается постоянным поперек полной толщины интерфейсного подслоя; 5(х - xf) - дельта-функция Дирака с Xf, являющейся мгновенным местоположением поверхности раздела. Дельта-функция появляется в вышеупомянутом выражении, потому, что оно было получено из
Vx = J§(х - Xf)nds.
Поверхностное натяжение о отражает избыток капиллярной энергии, сконцентрированной на поверхности раздела на единицу интерфейсной площади, вызванную изменением x через интерфейсный подслой.
divhI-gradh5p = -divhv или в операторном виде Ay = f
1
f « -| Vxl2 K Kn,
2
Монотонные численные схемы для уравнения переноса. Рассмотрим следующую схему [11]
[(У +1- У) /т] + а"у1 = 0, п = 0,1,... (7)
с нелинейным коэффициентом
(8)
_n n
a = a
n
1 + hyxx 2 n
v 2 Ух J
Нелинейная разностная схема (8) будет монотонной, если
ап > 0; (9)
тахаПТ ^ 1, п = 0, 1, .... (10)
1 ' к
Из (9), (10) следует, что нарушение монотонности наблюдается в окрестности экстремумов разностного решения.
Представим (8) в более удобном виде
~п п п /1
а = а х (11)
с множителем х
/ п п
п
X
1 + h УххУ.
2 I n\2 т 21 n
У^ + Yh Ух
(12)
и параметром Y- Дополнительные (регуляризирующие) слагаемые в %n из (12) имеют третий порядок по h.
Таким образом, регуляризованная схема (7), (11), (12) аппроксимирует одномерное уравнение адвекции в не дивергентной форме со вторым порядком по пространству.
Сформулируем условия монотонности схемы (7), (11), (12). Множитель a" (условие
(9)) будет неотрицательным при y = 0,25. Шаги сетки должны удовлетворять условию
(10), которое с учетом (11), (12) эквивалентно неравенству
maxanT( 1 + д1-) " n = 1,.... (13)
Таким образом, для регуляризованной схемы, (11), (12) при минимально допустимом
значении параметра a" максимально допустимый шаг по времени, согласно (13), уменьшается вдвое по сравнению со схемой с направленными разностями.
Предлагаемая вычислительная методика характеризуется высоким быстродействием, которое при решении CFD проблем равняется 5 х 10-6 с на узел на один шаг по времени. Что позволяет проводить вычисления на персональном компьютере с частотой 1 ГГц и RAM ~1 Гбайт.
Матрица валидации. Тест на давление внутри капли. В табл. 1 даны значения давления в пузырьке в зависимости от сетки в двумерном случае, в табл. 2 - значения давления в пузырьке в зависимости от сетки в трехмерном случае. Далее приведена матрица валидации по [5], для тестирования были отобраны 1, 2, 11 и 17 тесты.
Тест RISE OF A SPHERICA
Для дальнейшего прочтения статьи необходимо приобрести полный текст. Статьи высылаются в формате PDF на указанную при оплате почту. Время доставки составляет менее 10 минут. Стоимость одной статьи — 150 рублей.