CC BY
19
УДК 519.633, 532.529.6
ПРИМЕНЕНИЕ МЕТОДА ФУНКЦИИ УРОВНЯ В ЗАДАЧАХ ГИДРОДИНАМИКИ ДВУХФАЗНЫХ СРЕД
ТОНКОВ Л.Е.
Институт прикладной механики УрО РАН, 426067, г. Ижевск, ул. Т.Барамзиной, 34
АННОТАЦИЯ. Выполнено численное моделирование взаимодействия капель вязкой несжимаемой жидкости в среде идеального газа на основе метода функции уровня для системы уравнений Навье-Стокса. Рассмотрены процессы изменения поверхности раздела сред при дроблении и слиянии капель. Решения получены для фаз
с отношением плотностей 103 . Проведено сравнение с расчетными и экспериментальными данными других авторов.
КЛЮЧЕВЫЕ СЛОВА: численное моделирование, метод функции уровня, вязкая жидкость, свободная поверхность, уравнения Навье-Стокса.
ВВЕДЕНИЕ
Задачи, связанные с описанием двухфазных течений «газ» - «жидкость», во многих практически важных случаях можно свести к рассмотрению системы с несмешивающимися фазами. Моделирование такого рода течений приводит к необходимости решения системы уравнений Навье - Стокса для каждой фазы с учетом граничных условий на поверхности раздела и одновременным определением положения этой поверхности. Методы решения задач со свободными границами можно разделить на две группы. Значительная часть разработанных численных алгоритмов основана на процедуре явного выделения межфазных границ, когда расчетная сетка в процессе решения задачи постоянно перестраивается в соответствии с изменениями положения и формы поверхности раздела. Однако, описание слияния и разделения подобластей, занятых одной фазой (например, дробление струи на отдельные капли), как правило, приводит к существенному усложнению алгоритма.
Значительно более удобным может оказаться применение численных методов, в которых отсутствует необходимость явного выделения границы раздела путем перестроения расчетной сетки, то есть, методов сквозного счета. При этом полагается, что движение каждой фазы подчиняется одним и тем же уравнениям гидродинамики, но параметры среды (давление, плотность, вязкость) претерпевают разрыв на границе раздела фаз.
МАТЕМАТИЧЕСКАЯ МОДЕЛЬ
Пусть область О: О е Я заполнена двухфазной средой и О = и Г и О2 . Для определенности будем считать, что в подобласти О1 находится более плотная среда (жидкость), а в О2 - менее плотная (газ); Г является границей раздела фаз.
Для определения Г введем в рассмотрение скалярную функцию уровня x), (р\< 1 вид которой может быть в известном смысле произвольным, а сечение р(1 x) = 0 совпадает с границей раздела фаз [1]. Дополнительно к основной системе уравнений гидродинамики решается уравнение переноса для функции уровня р. Таким образом, формально в любой момент времени мы можем точно определить положение поверхности раздела. Будем далее рассматривать систему из двух несмешивающихся несжимаемых вязких сред, движение каждой из которых описывается уравнениями Навье-Стокса и уравнением неразрывности [2]
^^ + V • (рШ)г- = -Ур + V ■ тг- + д g, (1)
V- Ui= 0, (2)
где индекс / = {1,2} - номер среды, р - плотность, U - вектор скорости, р - давление, тi = / (VU + VU ) - тензор вязких напряжений, / - динамическая вязкость, g - вектор ускорения свободного падения.
На границе раздела сред Г выполняются условия динамического равновесия
(Т1 - = (Р1 - Р2 + оК)е, (3)
Ul = U2 , (4)
где e - единичный нормальный вектор, внешний по отношению к ^, К - кривизна поверхности Г , о - коэффициент поверхностного натяжения.
Геометрические характеристики поверхности раздела фаз Г однозначно определяются по функции уровня [1, с.131]
„ = к = vv?
V?? V? Г
Силу поверхностного натяжения, действующую в узком переходном слое, который в пределе является бесконечно тонким и совпадает с Г , можно свести к объемной [3]:
F = -о^дё?р), (5)
где ё - дельта-функция.
Обозначим через Н(?) функцию Хевисайда, такую, что Н(?) = 1 для ?> 0. Вводя в рассмотрение функции плотности р = р2 + (Р1 - Р2)Н(?), давления р = Р2 + (р - р2)Н(?) и вязкости / =/2 + (/1 - /2)Н(?), перепишем систему (1),(2) с учетом (3 - 5) следующим образом
дри + V-(рШ) = + V-т + рg + F, (6)
дt
V-U = 0. (7)
Система (6), (7) дополняется уравнением переноса функции уровня
д? = (8)
дt
ЧИСЛЕННЫЙ МЕТОД РЕШЕНИЯ
Для численного интегрирования системы (7 - 8) будем применять метод конечных объемов [4], который по определению обладает свойством консервативности и, кроме того, позволяет естественным образом использовать произвольные неструктурированные расчетные сетки, что делает его универсальным и удобным для расчета течений в областях со сложной геометрией.
Разделим область Q на непересекающиеся ячейки У : @ = и ¿У' и введем сеточные
функции /Р = (^,хг-), где ^ = Аt, п = 0,1,..., а хг - геометрический центр У . Для краткости под обозначением У будем понимать также значение объема г -й ячейки. Проинтегрируем (6) по У, обозначим через Н результирующий вектор объемных сил и с учетом теоремы Остроградского - Гаусса получим
— Г рШУ + ( (при)Ш? - ( = -( npds + ( + Г ШУ . (9)
Для вычисления производной по времени воспользуемся методом Эйлера первого порядка
^ ричУ . ^ - (ри)п У. (10)
дt1И Аt 1 '
Обозначим через Ni множество номеров ячеек сетки, являющихся соседними для Vi ,
то есть имеющих с ней общие грани. Вычисление диссипативных членов уравнения (10) в предположении ортогональности сетки сведется к простому суммированию
_ и - и
f n^VUds - £ И{ j-jS,, (11)
jeN. xj xi
где Sj = Sji , - площадь общей грани между ячейками i и j .
В случае неортогональной сетки для вычисления лапласиана (12) применяется какой-либо способ коррекции градиента, обеспечивающий второй порядок аппроксимации по пространству.
Наиболее существенным в рамках решения данной задачи является вычисление конвективных членов, так как именно конвективный перенос физических величин определяет общую картину течения, а наличие разрывов накладывает дополнительные ограничения на применяемые методы. Запишем конвективные слагаемые уравнения (10) в следующем общем виде
f (npUVdS - £ QVij. (12)
j^N.
Данное соотношение определяет конвективный перенос потоком Q некоторой физической величины ^ через грань между ячейками i и j . Для того, чтобы воспользоваться выражением (12) необходимо задать какой-либо способ вычисления (реконструкции) потоков Qj по известным значениям Qi и Qj . В данной работе для этого применена
численная схема MUSCL с ограничителем поток Ван-Лира [5]. Полностью аналогично решается и уравнение конвективного переноса функции уровня (8).
Подставим теперь (10 - 12) в (9), полагая, что интеграл от объемных сил вычисляется непосредственно по определению. Приводя далее подобные слагаемые, получим линейную систему алгебраических уравнений, которую будем решать методом сопряженных градиентов с предобуславливанием на основе неполной факторизации Холесского. Найденное таким образом поле скорости далее корректируется при помощи алгоритма PISO [6] для точного удовлетворения уравнению неразрывности (7).
В результате накопления в процессе решения вычислительных ошибок определенного вида (численной диффузии) функция уровня изменяет свой вид, что неизбежно приводит к погрешностям в определении величины р, а значит и к искажению геометрических характеристик поверхности раздела. Для управления влиянием численной диффузии возможно применение искусственного сглаживания разрывных функций на границе раздела фаз совместно с выполнением процедуры реинициализации [1, с.64]. Аналогично [7], функция Хевисайда H (р) заменялась приближенной.
РЕЗУЛЬТАТЫ
Все результаты получены методом конечных объемов на равномерной гексаэдральной сетке. Следует отметить, что описанный метод может быть применен без каких-либо изменений для любой неструктурированной сетки с произвольного вида ячейками при решении как двумерных, так и трехмерных задач.
Сравнение численного решения с результатами натурного эксперимента о соударении сферической капли и плоской преграды, выполненное в [8] показало хорошее совпадение рассчитанной картины изменения поверхности капли с данными скоростной киносъемки.
Рассмотрим в качестве первого примера задачу о неосевом соударении двух одинаковых капель вязкой жидкости [7]. На рис. 1 показаны основные фазы взаимодействия. Сначала происходит слияние капель, двигающихся навстречу с равными по модулю скоростями, в одну, которая, затем в силу закона сохранения импульса приобретает
вращательное движение вокруг центра масс = 0,06 с). Можно заметить, что центр масс рассматриваемой системы, являющийся также центром симметрии, покоится относительно наблюдателя (суммарный импульс равен нулю). Затем под влиянием сил инерции и поверхностного натяжения образуется характерный «перешеек» = 0,34 с), который в дальнейшем отделяется от основной массы жидкости = 0,44 с).
а) б)
а - расчет по представленной модели, б - результаты численного моделирования [7] Рис. 1. Динамика неосевого соударения двух капель
Различия представленных результатов с расчетом [7] в данном примере во многом объясняются достаточно грубой сеткой, выбранной из соображений уменьшения времени расчета. Размер к грани элемента равномерной кубической решетки задавался как h = К./30, где R — первоначальный радиус капли, что, по-видимому, является максимально допустимой величиной, обеспечивающей качественное согласование результатов моделирования и физической картины процесса.
Далее аналогично [9] рассмотрим задачу о соударении пяти капель вязкой несжимаемой жидкости отсутствии силы тяжести ^ = 0). В начальный момент времени капли уравновешенны внутренним давлением. Центральная капля покоится, а окружающие движутся по направлению к ее центру с равными по модулю скоростями. Отношение плотностей сред равно 1000, отношение вязкостей - 100. Число Рейнольдса Re = 100, число Вебера We = 100. Шаг равномерной расчетной сетки, как и в предыдущем примере к = Я/30.
На рис. 2 показаны сечения функции уровня ?(/,х) = 0, то есть конфигурация границы раздела жидкость-газ в различные моменты времени.
Рис. 2. Динамика плоскосимметричного соударения пяти одинаковых капель
Соударяясь, капли образуют единую систему (t = 2,0... 4,0) сложной формы. За счет кинетической энергии движущихся капель неподвижная ранее масса жидкости находящаяся в центре расчетной области начинает ускоряться в направлении, перпендикулярном плоскости движения (t = 6,0. 12,0). При этом за счет действия сил поверхностного натяжения и вязкости граница раздела становится боле гладкой, а ускорение, уменьшаясь, меняет свой знак на противоположный. Затем, образовавшаяся капля длительное время колеблется, сохраняя симметрию относительно центра расчетной области (t = 40,0). Через достаточно продолжительный промежуток времени силы поверхностного натяжения приводят форму границы раздела к сферической (на рисунке не показано).
ЗАКЛЮЧЕНИЕ
Метод функции уровня в сочетании с конечно-объемной численной схемой интегрирования системы уравнений гидродинамики позволяет моделировать сложные течения многофазных сред с изменяющейся топологией поверхностей раздела с достаточной для многих практических приложений с точностью без привлечения дополнительных процедур уточнения положения свободных границ. При этом обеспечивается соблюдение законов сохранения динамики сплошных сред, в том числе для течений, содержащих разрывы гидродинамический величин.
СПИСОК ЛИТЕРАТУРЫ
1. Osher S.J., Fedkiw R.P. Level Set methods and dynamic implict surfaces. Berlin, Springer, 2003. 273 p.
2. Лойцянский Л.Г. Механика жидкости и газа. М. : Наука, 1987. 840 с.
3. Brackbill J.U., Kothe D.B., Zemach C.A. А ^ntinuum method for modelling surface tension // J. Comp. Physics. 1992. V.100, №2. P.335-354.
4. Ильин В.П. Методы конечных разностей и конечных объемов для эллиптических уравнений. Новосибирск : Изд-во Ин-та математики СО РАН, 2002. 345 c.
5. Van Leer B. Towards the ultimate conservative difference scheme. II. Monotonicity and conservation combined in a second-order scheme // J. Comp. Physics. 1977. V. 32, №4. P.361-370.
6. Issa R.I. Solution of the implicitly discretized fluiduid flow equations by operator-splitting // J. Comp. Physics. 1986. V. 62, №1. P. 40-65.
7. Dai M., Schmidt D.P. Numerical simulation of head-on droplet collision: Effect of viscosity on maximum deformation // Phys. Fluids-2005. V. 17, №4. P. 041701.1-041701.4.
8. Тонков Л.Е. Численное моделирование динамики капли вязкой жидкости методом функции уровня // Вестник Удмуртского университета. Математика. Механика. Компьютерные науки. 2010. Т.3. С. 134-140.
9. Liang R., Satofuka N. Numerical study on three-dimensional unsteady motion of droplets using level set method // Proc. 4th ECCOMAS conference «Computational Fluid Dynamics-98». Greece, Athens : John Wiley & Sons, 1998. V.2. P.421-425.
LEVEL SET METHOD FOR TWO PHASE HYDRODYNAMIC PROBLEMS
Tonkov L.E.
Institute of Applied Mechanics, Ural Branch of the Russian Academy of Sciences, Izhevsk, Russia
SUMMARY. In this paper, the numerical simulation of viscous drop dynamics was studied by level set method for the incompressible Navie-Stokes equations. Solution procedure employs finite volume method on unstructured hexahedral grid elements. Some numerical results are presented and compared with other simulations.
KEYWORDS: numerical simulations, level set method, free surface tension, Navier-Stokes equations.
Тонков Леонид Евгеньевич, кандидат физико-математических наук, старший научный сотрудник ИПМ УрО РАН, тел. (3412) 20-35-14, e-mail: tnk@udman.ru
