CC BY
318
УДК 532.582.8
Ю.А. Крюков1
МОДЕЛИРОВАНИЕ ДВИЖЕНИЯ СФЕРИЧЕСКОЙ КАПЛИ В ЛАМИНАРНОМ ПОГРАНИЧНОМ СЛОЕ БЛАЗИУСА С ПОМОЩЬЮ ПАКЕТА ANSYS FLUENT
В статье описывается решение задачи о движении абсолютно твердой сферической капли малых размеров в ламинарном пограничном слое плоской полубесконечной пластины (Блазиуса) без учета ее вращения в пакете Ansys Fluent. Дается описание построения расчетной области задачи, создание структуры сеточных моделей, ставятся граничные и начальные условия на границах области и на капле, приводятся последовательность и ход решения.
Рассмотрены примеры решения задачи при различных начальных условиях капли в сравнении с методикой, основанной на интегрировании уравнений движения капли с учетом заранее определенных активно действующих на нее сил. Получено удовлетворительное совпадение результатов. Сделан вывод о целесообразности сочетания различных подходов для успешного решения подобных задач.
Ключевые слова: ламинарный пограничный слой, сферическая капля, задача Блазиуса, сеточный генератор ICEM CFD, динамически перестраивающиеся сетки, Ansys Fluent, решатель 6 DOF, уравнения движения Ньютона, активно действующие силы, поперечная сила, сила Сэфмана.
Введение
Для изучения движения тела в потоке жидкости или газа широко распространен подход, основанный на применении уравнений движения Ньютона. Существует два способа реализации этого подхода. Первый заключается в том, что исследователь выбирает систему сил, активно действующих на тело, записывает уравнение движения тела и интегрирует их. Основными вопросами при использовании этого подхода являются выбор системы активных сил и влияние тела на несущую фазу [1; 2]. Второй способ избавляет исследователя от решения указанных проблем, так как при этом рассматривается совместное движение тела и потока. Целью настоящей статьи является описание применения второго способа на примере движения сферической капли в ламинарном пограничном слое (ЛПС) плоской пластины с использованием пакета Ansys Fluent и сравнение полученных результатов с результатами, полученными при использовании первого способа.
!© Крюков Ю.А., 2015
Крюков Юрий Александрович (yurakryukov1985@mail.ru), кафедра математического моделирования в механике, Самарский государственный университет, 443011, Российская Федерация, г. Самара, ул. Акад. Павлова, 1.
1. Физическая постановка задачи
Рассмотрим движение сферической капли в ЛПС Блазиуса плоской полубесконечной пластины (рис. 1).
Рис. 1. Схема движения капли:
1 — пластина, 3 — ЛПС, Oxy — декартова система координат, (x0,y0) — координаты начальной точки движения капли в ЛПС, (xi,0) — координата конечной точки ее движения, (x|,0) — координата конечной точки капли без учета ЛПС, а — угол наклона пластины к горизонту, иж — скорость набегающего потока, g — вектор ускорения свободного падения
Выберем декартову систему координат и поместим ее на краю пластины. Пусть пластина, наклоненная под углом а к горизонту, обтекается невозмущенным потоком воздуха с постоянной скоростью иж, вектор которой совпадает с плоскостью пластины. Сферическая капля, свободно падая под действием сил тяжести, со скоростью Уж попадает в область ЛПС (точка с координатой (x0,y0) на рис. 1) и движется далее, пока не достигнет пластины (точка (xi,0)). Будем считать, что капля не деформируется и не вращается. Определим движение капли в ЛПС. Для решения поставленной задачи применим пакет вычислительной гидромеханики Ansys Fluent: поместим каплю в область ЛПС и зафиксируем ее положение. Сформируем поле течения, для чего найдем стационарное решение задачи. Затем перейдем к нестационарному решению, задавая начальные условия для капли. Будем следить за перемещением капли, используя встроенный решатель 6 DOF. Этот решатель определяет перемещение тела с учетом приложенных на него массовых и объемных сил, которые получает численным интегрированием давления и касательного напряжения по поверхности движущегося тела. Перемещение капли в границах расчетной области будет возможным благодаря использованию подхода динамически перестраивающейся сетки, который позволяет передвигать границы ячеек и деформировать сетку с учетом движения капли [3; 4].
2. Построение расчетной области
Для расчета движения капли в Ansys Fluent нет необходимости моделировать весь процесс, начиная от свободного падения и заканчивая касанием каплей пластины (это может значительно увеличить время расчета). Скорость Уж, с которой
X
капля попадает на внешнюю границу ЛПС, может быть разной: меньше или равной скорости витания Vmax, если капля свободно падает, или больше Vmax, если капля получила дополнительное ускорение (например, при соударении с другой каплей). Скорость витания определим из баланса сил тяжести и сопротивления, действующих на нее при свободном падении в невозмущенном воздухе [5]. Для этого в пакете Ansys Fluent выполним последовательность расчетов по обтеканию сферической капли (например, с учетом методики, описанной в [6]). Расчеты будем продолжать до тех пор, пока сила сопротивления не сравняется с силой тяжести. Та скорость набегающего потока, при которой сила сопротивления, действующая на каплю, уравновешивает силу тяжести, будет скоростью витания.
Расчетную область и сетку будем строить с использованием генератора сеточных моделей ICEM CFD [4]. Всю расчетную область представим в виде криволинейной трапеции той области ЛПС, в границах которой будет двигаться капля (рис. 2).
Рис. 2. Расчетная область:
ивь, ^еь — соответственно профиль продольной и поперечной скорости ЛПС, V — кинематическая вязкость
Левую и правую границы области ограничим продольными координатами хо — —Дх и хо+Д2 соответственно. Очевидно, что Д1 и Д2 нужно выбирать так, чтобы капля во все время движения оставалась в расчетной области, то есть Д1 > 0 и Д2 > (х\ — хо). Значение х\ можно определить, решив в подвижных координатах задачу: пусть капля с начальными координатой (хо,уо) и компонентами горизонтальной иж и вертикальной —Утах скорости падает на пластину, наклоненную под углом а к горизонту. Примем за единицу время полета капли до пластины, а влияние ЛПС на каплю учитывать не будем. Тогда справедливо соотношение х—Х0 — , из которого выразим координату столкновения капли с пла-
У0 Vmax cos a
стиной x\ = Vmm y0 + x0. Верхнюю границу определим по формуле Блази-
уса для толщины пограничного слоя 5 = 6Л/[7]. Граничные условия опреде-
у u<Xj
лим, используя UDF-макрос [8]. На левой границе зададим известные профили продольной и поперечной скоростей Блазиуса при координате хо — Д1, на верхней границе - постоянную по величине продольную скорость, равную иж, а также профиль поперечной скорости vbl(5) = 0, 857y/V3X~ , определенный по известной формуле Блазиуса на внешней границе ЛПС [9]. На левой и верхней границах поставим условие типа "inlet". На поверхности пластины и для капли — условия прилипания и непротекания потока (типа "wall"), на правой границе — условие
типа "pressure outlet", определяющее статическое давление на выходе потока из области. На рис. 2 возле границ расчетной области указан тип соответствующих граничных условий и профили граничных скоростей.
Рассматриваемая задача является трехмерной, плотности фаз сильно различаются, то есть рз ^ р2 (индекс 2 относится к жидкости). Согласно предварительным расчетам, для капли можно допустить наличие только двух степеней свободы (капля перемещается в направлении х и у). Поэтому в направлении z (на рис. 2 ось z перпендикулярна плоскости оху) длину расчетной области можно оценить порядком нескольких диаметров капли. Каплю диаметром D2 поместим в точку с координатой (xo, yo = 5^/Ц^0). На боковых поверхностях расчетной области (параллельных плоскости оху) зададим граничное условие типа "symmetry" (симметрия).
3. Построение расчетной сетки
Построим в расчетной области сетку, состоящую из двух зон. Одна из них примыкает непосредственно к капле (зона "inner"), в процессе расчета не деформируется и перемещается вместе с ней. Она состоит из гексаэдральных ячеек. На рис. 3 представлен фрагмент расчетной сетки в плоскости оху миделевого сечения капли.
Зона "inner" построена в виде полого шара, внутренней поверхностью она описывает сферическую поверхность капли, а на внешней границе плавно переходит в зону "outer". Зона "outer", состоящая из тетраэдральных ячеек, заполняет оставшуюся расчетную область и имеет большие объемы ячеек, чем зона "inner". Это необходимо для лучшего разрешения уравнений в области пограничного слоя капли [6]. При расчете ячейки зоны "outer" перестраиваются, деформируются и разбиваются с учетом движения зоны "inner". Окончанием расчета является время, когда зона "inner" достигнет пластины.
4. Пример
Приведем пример решения задачи с помощью описанной методики для капли диаметром D2 = 10-4 м (состав — вода). Пусть окружающая среда — воздух при н. у. Поместим каплю в ЛПС (иж = 5 м/с ) горизонтальной пластины в точку с абсциссой xo = 1 ми ординатой yo =8, 63 • 10-3 м . По приведенным данным определим число Рейнольдса для капли Кв2 = 33, 6 (за характерную скорость взята иж). Тогда, согласно [6], ламинарную модель течения также можно использовать для решения задачи об обтекании капли. Получим стационарное решение задачи, результатом которой будет сформировавшееся поле скоростей в рассматриваемой области. Далее перейдем к нестационарной задаче, для чего назначим начальные условия для капли: начальная продольная скорость капли xo = 5 м/с, поперечная yo = Vmax = _0, 2495 м/с .
Для перестраивающийся сеточной зоны "outer" используются методы сглаживания "smoothing" с настройками по умолчанию и перестроения "remeshing" с опцией "local cell" и настройками согласно рекомендациям [3]. Общее число ячеек во всей области на начальный момент времени n = 5 • 106 шт., а все время расчета занимает порядка 2-3 суток с использованием 128 ядер на кластере.
Интерес представляет соотношение сил, действующих на каплю и определяющих ее движение. Введем в рассмотрение коэффициенты объемных продольной Cx и поперечной Cy сил, а также коэффициент силы тяжести Cg . Все коэффициенты отнесены к площади миделевого сечения капли и скоростному напору, где в качестве характерной скорости взята иж. Для оценки достоверности полученных результатов на некоторых графиках, приведенных ниже, сплошной линией приведены кривые, полученные с помощью второго способа (Ansys Fluent), а пунктирной — с помощью первого способа (Maple).
Результаты и выводы
На рис. 4 представлено изменение по времени коэффициента продольной силы, действующей на каплю, и ее продольной скорости х.
Рис. 4. Коэффициент продольной силы и продольная скорость капли
На рис. 4 по графику слева видно, что сила является отрицательной (оказывает сопротивление движению капли), и с приближением капли к пластине ее
величина монотонно увеличивается до экстремальных значений. Скорость капли при этом, очевидно, падает (на рис. 4 график справа).
Изменение по времени коэффициентов сил и поперечной скорости у капли приведены на рис. 5.
Рис. 5. Коэффициенты поперечной силы, силы тяжести и поперечная скорость капли
По рисунку видно, что Су является положительным, и соответствующая сила стремится оттолкнуть каплю от пластины. По величине Су превосходит Сд (график слева на рис. 5), поэтому поперечная скорость капли падает (график справа на рис. 5).
Изменим начальные условия для капли. Пусть она попадает в ЛПС со скоростью меньше равновесной у о = Уж = —0,14 м/с. На рис. 6 приведены коэффициенты сил и поперечная скорость капли у.
Рис. 6. Коэффициенты сил и поперечная скорость капли при неравновесной начальной скорости
По графикам заметно, что баланс сил, действующих на каплю в поперечном направлении, меняется. Так, при Ь < 0,23 с величина Сд превосходит Су (график слева на рис. 6), и капля продолжает набирать скорость (график справа на рис. 6). При Ь > 0,23 с соотношение сил изменяется на обратное, и капля тормозится.
На рис. 4-6 (графики слева) асимптотическое поведение кривых в окрестности начала отсчета по времени объясняется особенностью решения задачи в Ansys Fluent при переходе от стационарного решения к нестационарному.
В пакете Ansys Fluent объемная сила, действующая на каплю со стороны потока, получается по результатам интегрирования распределений давления и касательного напряжения по поверхности капли. Таким образом, эта сила является суммарной, и расчленить ее на сумму сил разной природы в этом пакете затруднительно. Гораздо проще продемонстрировать соотношение сил, используя Maple. На рис. 5 и 6 (графики слева) приведены коэффициенты поперечной силы сопротивления Cd и силы Сэфмана Csaf сдвиговой природы [10-14], суммарно равные коэффициенту объемной силы Cy. По рисункам видно, что величины этих сил одного порядка, поэтому при моделировании движения капли указанные силы стоит учитывать.
Сравнение результатов решения, полученных с помощью Ansys Fluent и Maple, которое можно сделать по рис. 4-6, показывает их хорошее совпадение. Таким образом, представленную задачу можно решать, используя как первый, так и второй способы. Первый способ отличается быстрым получением результатов (расчет систем обыкновенных дифференциальных уравнений занимает время порядка минуты), но требует проверки правильного выбора системы активно действующих сил. Второй способ много затратнее первого по времени и ресурсам, однако избавляет исследователя от выбора системы сил и дает более точные результаты, так как еще учитывает взаимное влияние несущей среды и тела. Поэтому для успешного решения задачи о движении тела целесообразным видится подход, когда применяются оба способа: с помощью результатов, полученных в пакете Fluent, можно проверять систему активных сил, используемых в способе Maple, тем самым отрабатывая способ Maple для серийных расчетов.
Литература
[1] Рыбдылова О.Д. Поперечная миграция частиц в пограничном слое на плоской пластине // Наука и образование: эл. научно-техн. изд. 2012. № 2. 77-30569/318597. URL:http://technomag.edu.ru/doc/318597.html
[2] Клюев Н.И., Мингулов Х.И., Бурмистров Н.А. Движение капли в градиентном потоке // Вестник СамГУ. Естественнонаучная серия. 2012. № 3/2(94). С. 24-28.
[3] ANSYS FLUENT 12.0 Theory Guide, 2009.
URL: http://orange.engr.ucdavis.edU/Documentation12.0/120/FLUENT/flth.pdf (дата обращения: 19.01.2015)
[4] ANSYS ICEM CFD User Manual, 2012. URL:http://ptechgroup.ir/Virtual Education/Mechanics/ICEM-CFD/PTEC-ICEM CFD 14.5 User Guide.pdf (дата обращения: 19.01.2015)
[5] Уоллис Г., Одномерные двухфазные течения. М.: Мир, 1972.
[6] Jones D.A., Clarke D.B. Simulation of Flow Past a Sphere using the Fluent Code, 2008. URL:http://www.dsto.defence.gov.au/corporate/reports/DSTO-TR-2232.pdf (accesed: 22.10.2014)
[7] Шлихтинг Г. Теория пограничного слоя. М.: Наука, 1969. 744 c.
[8] ANSYS FLUENT 12.0 UDF MANUAL, 2009.
URL: http://orange.engr.ucdavis.edu/Documentation12.0/120/FLUENT/fludf.pdf (дата обращения: 19.01.2015)
[9] Лойцянский Л.Г. Ламинарный пограничный слой. М.: Физматгиз, 1962. 480 с.
[10] Saffman P.G. The lift on a small sphere in a slow shear flow // J. Fluid Mech. 1965. V. 22. P. 385-400.
[11] McLaughlin J.B. Inertial migration of a small sphere in linear shear flows //J. Fluid Mech. 1991. V. 224. P. 261-274.
[12] McLaughlin J.B. The lift on a small sphere in wall-bounded linear shear flows // J. Fluid Mech. 1993. V. 246. P. 249-265.
[13] Mei R. An approximate expression for the shear lift force at finite Reynolds number // Intern. J. Multiphase Flow. 1992. V. 18. N 1. P. 145-147.
[14] Multiphase flows with droplets and particles / C.T. Crowe [et al.]. CRC Press. 2012. 487 p.
References
[1] Rybdylova O.D. Transverse migration of particles in the boundary layer on a flat plate]. Nauka i Obrazovanie: el.nauchno-technich. izdanie [Science and Education: Scientific and Technological Publication], 2012, no. 2, 77-30569/318597. Retrieved from: http://technomag.edu.ru/doc/318597.html (accesed 19.01.2015) [in Russian].
[2] Klyuev N.I., Mingulov Kh.I., Burmistrov N.A. The movenment of a drop in the gradient flow. Vestnik Samarskogo gosudarstvennogo universiteta. Estestvenno-nauchnaya seriya [Vestnik of Samara State University. Natural Science Series], 2012, no. 3/2(94), pp. 24-28 [in Russian].
[3] ANSYS FLUENT 12.0 Theory Guide, 2009. Retrieved from: http://orange.engr.ucdavis.edu/Documentation12.0/120/FLUENT/flth.pdf (accesed 19.01.2015)
[4] ANSYS ICEM CFD User Manual, 2012. Retrieved from: http://ptechgroup.ir/Virtual Education/Mechanics/ICEM-CFD/PTEC-ICEM CFD 14.5 User Guide.pdf (accesed 19.01.2015)
[5] Uollis G. One-dimensional two-phase flows. M.,Mir, 1972, 440 p. [in Russian].
[6] Jones D.A. and Clarke D.B. Simulation of Flow Past a Sphere using the Fluent Code, 2008. Retrieved from: http://www.dsto.defence.gov.au/corporate/reports/DSTO-TR-2232.pdf (accesed
22.10.2014)
[7] Shlihting G. Boundary layer theory. M., Nauka, 1969, 744 p. [in Russian].
[8] ANSYS FLUENT 12.0 UDF MANUAL, 2009. Retrieved from: http://orange.engr.ucdavis.edu/Documentation12.0/120/FLUENT/fludf.pdf (accesed
19.01.2015)
[9] Loitsiaynsky L.G. Laminar boundary layer. M., Fizsmatgiz, 1962, 480 p. [in Russian].
[10] Saffman P.G. The lift on a small sphere in a slow shear flow. J. Fluid Mech., 1965, Vol. 22, pp. 385-400.
[11] McLaughlin J.B. Inertial migration of a small sphere in linear shear flows. J. Fluid Mech., 1991, Vol. 224, pp. 261-274.
[12] McLaughlin J.B. The lift on a small sphere in wall-bounded linear shear flows. J. Fluid Mech., 1993, Vol. 246, pp. 249-265.
[13] Mei R. An approximate expression for the shear lift force at finite Reynolds number. Intern. J. Multiphase Flow, 1992, Vol. 18, no. 1, pp. 145-147.
[14] Crowe C.T., Schwarzkopf J.D., Sommerfeld M., Tsuji Y. Multiphase flows with droplets and particles. CRC Press, 2012, 487 p.
Yu.A. Kryukov2
MODELING THE MOTION OF A SPHERICAL DROP INTO A LAMINAR BOUNDARY LAYER USING ANSYS FLUENT
The article provides a formulation to the problem of spherical drop motion into laminar boundary layer over a flat semi-infinite plate (Blasius problem). The paper presents the solution in the package Ansys Fluent without taking into account deformation and rotation of drops. The article describes building of geometric area, construction of the grid models, boundary and initial conditions statement, methods of solution, course of computation. The article describes the examples of problem solving for various initial conditions of a drop, compares the results obtained using different methods. There is good agreement between the results. A combination of different methods for solving similar problems is the guarantee of a successful solution.
Key words: laminar boundary layer, spherical drop, Blasius task, ICEM CFD, dynamic meshes, Ansys Fluent, 6 DOF Solver, Newton's equation of motion, active forces, transverse force, Saffman lift force.
Статья поступила в редакцию 19////2014. The article received 19////2014.
2 Kryukov Yuriy Alexandrovich (yurakryukov1985@mail.ru), Department of Mathematical Modelling in Mechanics, Samara State University, 1, Acad. Pavlov Street, Samara, 443011, Russian Federation.
