Расчет поля скорости при течении вязкоупругой жидкости в аппарате с мешалкой Текст научной статьи по специальности «Физика»

УДК 532.517.2

РАСЧЕТ ПОЛЯ СКОРОСТИ ПРИ ТЕЧЕНИИ ВЯЗКОУПРУГОЙ ЖИДКОСТИ В АППАРАТЕ С МЕШАЛКОЙ

© 2014 Н А. Газизуллин

Казанский национальный исследовательский технологический университет

Поступила в редакцию 01.03.2014

Методом контрольных объемов выполнено численное моделирование ламинарного течения вязко-упругой жидкости в аппарате с механической мешалкой. Результаты расчетов представлены в виде линий тока вторичной циркуляции жидкости.

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

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

В данной работе проведено численное моделирование ламинарного течения несжимаемой неньютоновской жидкости, описываемой реологической моделью вязкоупругой Олдройд-Б жидкости, в аппарате с турбинной мешалкой с вертикальными прямыми лопатками полной длины, расположенными радиально (рис. 1). Основное преимущество турбинных мешалок состоит в возможности их применения в широком диапазоне изменения вязкости и плотности перемешиваемых сред.

Исходными уравнениями, описывающими течение жидкости в аппарате, будут уравнения движения и неразрывности в виде [4]

д ^ —(pv) + V ■ pvv = F + V ■ (- pI + а),

dt (1) V- = 0, (2)

где p - плотность жидкости, p - давление, о -тензор-девиатор напряжений, I - единичный тензор.

Газизуллин Назым Абдуллович, кандидат технических наук, доцент кафедры высшей математики.E-mail: gnazym @gmail. com

Рис. 1. Турбинная мешалка

Реологическая модель вязкоупругой Олдройд-Б жидкости записывается в виде соотношения [5]

а + \ а = 2ц

D + D

V

У

(3)

где п - сдвиговая вязкость, - время релаксации, А,2 - время ретардации, а тензор скоростей деформаций определяется как

D=1И+(Уу )Т ].

Соотношение (3) содержит также верхнюю конвективную производную, определяемую в виде [5]

а = — + V • Уа - (Уу )Т ■ а - а ■ Уу ,

&

где верхний индекс Т обозначает операцию транспонирования.

V

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

¥ = р(ю2 г + 2юv),

¥Ф=-2Рюи , ¥г =-РЯ.

Для решения системы уравнений (1), (2) и (3) используем метод расщепления напряжений [6, 7], согласно которому тензор напряжений расщепляется на вязкоупругую составляющую и вязкую составляющую равенством

а = т + 2^2 Б,

(4)

где вязкоупругая составляющая тензора напряжений т удовлетворяет уравнению

т + ^т = 2^ Б.

(5)

В формуле (4) постоянные п и п2, представляющие собой соответственно неньютоновский и ньютоновский вклады в общую вязкость, определяются соотношениями

Л = Л1 + Л2, ^ =(1 -Р)л^1, Р=—

Л

С учетом (4) уравнение движения (1) примет вид

д_ дг

(рV) + V • (рVV - лVV) =

= ¥ -Vp + V • (т - 2рлБ).

(6)

Граничные условия на свободной поверхности включают в себя кинематическое и динамическое условия. Кинематическое условие

дИ дИ

— = -и--+ w,

дг дг

(7)

достаточно велик. Тогда динамическое условие может быть записано в виде

п

(- р1 + а) = -Роп • 1,

(8)

где ро - атмосферное давление, а вектор нормали к свободной поверхности жидкости

п-{-! >4

Граничные условия для составляющих скорости на твердых стенках заключаются в отсутствии относительного движения жидкости и твердой поверхности. Тогда на дне и боковой стенке аппарата

и = 0, V = —юг, w = 0, а на поверхности вала и мешалки соответственно

и = 0, V = 0, w = 0,

где ю - угловая скорость вращения вала и мешалки.

На оси вращения потока под мешалкой примем

дw

и = 0, V = 0, — = 0, дг

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

* иг * г *

г =—, г =—, ф =ф, Ь V

*

г =<

Ь

г - Н - Ь

И

при г < Н + Ь при г > Н + Ь

где И - высота свободной поверхности над мешалкой, свидетельствует о том, что скорость движения свободной поверхности в направлении нормали к ней совпадает с нормальной составляющей скорости движения жидкости [8]. Динамическое условие является условием равенства сил, действующих на свободную поверхность [9]. Будем предполагать, что силы поверхностного натяжения незначительны, так как радиус кривизны поверхности жидкости в аппарате

- = 1 * = Ь(Р - Р0 )

и' ли ' ли

* Ьа а =-

где Н - высота расположения мешалки над дном аппарата; Ь - высота лопасти мешалки, а в качестве характерной длины Ь и характерной скорости потока и выбраны соответственно диаметр мешалки ё и окружная скорость (в неподвижной системе координат) конца лопасти. Заметим, что свободной поверхности жидкости соответствует

* 1

при этом значение г =1.

V

г

После преобразования координат уравнение неразрывности (2) сохранит вид

V-V = 0.

При этом проекции вектора скорости в расчетной области будут определяться следующим образом:

и = уи *, V = уу*, Ж = м>* - и* г * ,

дг

где у=1 для нижней подобласти и у=Н* для верхней подобласти; И =И/ё, а и , у , w - компоненты вектора скорости в безразмерной физической области, определяемой преобразованием

* г * * г

г =-, ф =Ф, * =-

а а

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

д/ у— д хк

И1И2 И3 ^ п 2 3 Ли, Ф

И

1

д

' И1И2И3 дФ

— дX

И,2

д х

+ &

к у

(9)

полшага в соответствующих направлениях относительно основных точек, в которых вычисляются давление и компоненты тензора напряжений. Проинтегрируем уравнение (9) по контрольному объему и временному интервалу А/. В результате, с учетом уравнения неразрывности, получим соотношение, связывающее значения искомой функции Ф в узловой точке Р с ее значениями в центрах Е, Ж, Ы, £, Т, В соседних ячеек, в виде

арФ р = аЕ Ф е + а№ Ф№ + Ф N + +а Ф + а Ф + а? Ф + AV,

(10)

где £Р - узловое значение источникового члена; А^- объем ячейки.

Для расчета формы свободной поверхности жидкости используем динамическое условие (8) в проекции на нормаль

„ дИ

агг - Р + Р0 - 2аг^ +

+ (агг - Р + Ро )1д- I = 0

дИ

дг

2

дг

(11)

где приняты следующие обозначения: -1=1, И2=г , -3=1, х1=г , х2=ф , х3=г , и1=и , и2=у , и3=^ , к - индекс суммирования; Г=1; Л=пЯе; Яе=рпс!1/п - центробежное число Рейнольдса; п - число оборотов мешалки в единицу времени, а £Ф - член типа источника, который определяется соответственно искомой функции Ф.

Уравнение (5) в координатной форме также может быть представлено в виде обобщенного уравнения (9). При этом следует считать Г=0; Л=Же; Же=А1ю - число Вайссенберга. Численное моделирование течения проведем методом контрольных объемов [10], одним из основных достоинств которого является то, что он обладает консервативными свойствами, то есть обеспечивает интегральное выполнение законов сохранения для всей расчетной области. Поделим расчетную область на контрольные объемы (ячейки) таким образом, чтобы каждая узловая точка находилась в отдельной ячейке. При этом грани контрольных объемов расположим посередине между соседними узловыми точками. Размещение всех узловых функций в одних и тех же точках приводит к рассогласованию полей скорости и давления, поэтому выберем разнесенную шахматную сетку [11], в которой точки, где вычисляются компоненты скорости, смещены на

где огг, о22, о^ - компоненты тензора напряжений.

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

- (/ *, г* )= -* (/'

', г* +Аг *), И*(/*,Я*)= И*(/*,Я* -Аг*),

(12) (13)

где гв - радиус вала; Я - радиус аппарата. Перепишем условие (11) в безразмерном виде

аИ*

дг *

а гг - Р +(а гг - Р

(* * \ а*г - Р )

2а*

(14)

Решение уравнения (14) проводилось на каждой итерации методом Рунге-Кутта [12] при начальном условии И0 =И -А/ ,гв +Аг) с учетом условий (12), (13) при фиксированном значении производной дИ /дг , взятом с предыдущей итерации. При вычисленных на к-ой итерации значениях высоты свободной поверхности И* *(кк может не выполняться условие постоянства объема жидкости в аппарате (при этом достаточно учитывать объем жидкости над мешалкой). Таким образом, возникает необходимость введения некоторой поправки 5Ик к значениям И**(кк на

2

ж

каждом итерационном шаге. Эта поправка находилась из соотношения

V - V

R*2 - г;2)h\

(15)

h4 k ) = h**( k ) +8h(k ).

(16)

aP bpP = aE bpE + aw bpw + aN bpN + + as Sps + a Sp + aB SpB + bp

ue = u* +

* ** u = u +

w w

* ** v„ = v„ +

ГеАффА- (Spp-Spe ),

ae

rwA^Az* (Spw-SPP ),

(5pp -Spn ),

a

Ar * Az

w

* «

* ** ArAz /„ ,,4 v* = v* +-№ - Spp Л

5 5

* ** w* = w* +

a

s

r Аф Ar (spp -Spr ),

где ¥0 - объем жидкости над мешалкой с невозмущенной свободной поверхностью; V® - объем жидкости над мешалкой на к-ой итерации, который вычислялся на каждой итерации путем чис-

1 **(к)

ленного интегрирования по значениям к ' Л Следует отметить, что в соответствии с формулой (15) поправка 5кк не влияет на форму свободной поверхности, а лишь корректирует ее по высоте. Окончательно с учетом поправки скорректированные значения высоты свободной поверхности к**® могут быть найдены как

lat

* ** ! г.

wb = w* +

Г Аф Az / ч --(Spb -spp л

yab

и давления

p* = p** + aôp

Расчет поля течения проводился на основе алгоритма SIMPLE [10], в котором используется дискретизация уравнений по методу контрольных объемов на сетках с расположением узлов в шахматном порядке. В алгоритме SIMPLE должно задаваться начальное предполагаемое поле

давления. Это поле давления в соответствии с

**

рекомендациями [10] полагалось равным p =0 для того, чтобы погрешности округления при расчете градиентов давления не достигали слишком больших значений. Вначале по уравнениям (10) вычислялись компоненты тензора напряжений и предварительные значения компонент скорости u , v , w , не удовлетворяющие уравнению неразрывности. С учетом приближенного решения для скорости находилась поправка к давлению 5p из уравнения

(17)

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

Затем с учетом поправок 5р рассчитывались скорректированные значения компонент скорости на гранях е, w, п, s, г, Ь ячеек, удовлетворяющие уравнению неразрывности, по поправочным формулам

где a - параметр релаксации. В расчетах принималось a=0,8.

Форма свободной поверхности жидкости, соответствующая рассчитанному полю течения

определялась в конце каждой итерации, при

" 1 *(0)

этом в качестве начальных значений h,y ' принималось значение, соответствующее положению невозмущенной поверхности жидкости. Дискретные уравнения (10) и (17) решались методом прогонки [13]. В качестве критерия сходимости рассматривалась сумма модулей невязок по всем контрольным объемам при решении уравнений (10). Расчеты проводились по этому критерию с точностью до 10"6 на равномерной сетке. В расчетах принималось H0=D, d/D = 0,3, de /D = 0,05, b/d = 0,2, H/H0 = 0,4, где D - диаметр аппарата; H0 - высота невозмущенной поверхности жидкости над дном аппарата; dB - диаметр вала.

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

n

2.

3.

Рис. 2. Радиально-осевая циркуляция в аппарате при Яе = 150, Же = 5.

СПИСОК ЛИТЕРАТУРЫ:

Стренк, Ф. Перемешивание и аппараты с мешалками. - Л.: Химия, 1975. 384 с. Брагинский, Л.Н. Перемешивание в жидких средах / Л.Н. Брагинский, В.И. Бегачев, В.М. Барабаш. -Л.: Химия, 1984. 336 с.

Манусов, Е.Б. Расчет реакторов объемного типа / Е.Б. Манусов, Е.А. Буянов. - М.: Машиностроение, 1978. 112 с.

4. Астарита, Дж. Основы гидромеханики неньютоновских жидкостей / Дж. Астарита, Дж. Маруч-чи. - М.: Мир, 1978. 312 с.

5. Larson, R.G. The structure and rheology of complex fluids. - New York: Oxford University Press, 1999. 663 pp.

6. Xue, S.C. Three-dimensional numerical simulations of viscoelastic flows through planar contractions / S.C. Xue, N. Phan-Thien, R.I. Tanner // Journal of Non-Newtonian Fluid Mechanics. 1998. Vol. 74. P. 195245.

7. Phillips, T.N. Viscoelastic flow through a planar contraction using a semi-Lagrangian finite volume method / T.N. Phillips, A.J. Williams // Journal of Non-Newtonian Fluid Mechanics. 1999. Vol. 87. P. 215246.

8. Лаврентьев, МА. Проблемы гидродинамики и их математические модели / МА. Лаврентьев, Б.В. Шабат. - М.: Наука, 1977. 408 с.

9. Ландау, Л.Д. Гидродинамика / Л.Д. Ландау, Е.М. Лифшиц. - М.: Физматлит, 2003. 732 с.

10. Патанкар, С. Численные методы решения задач теплообмена и динамики жидкости. - М.: Энерго-атомиздат, 1984. - 152 с.

11. Harlow, F.N. Numerical calculation of time-dependent viscous incompressible flow of fluid with free surface / F.N. Harlow, J.E. Welch // Physics of Fluids. 1965. Vol. 8. P. 2182-2189.

12. Холл, Д. Современные численные методы решения обыкновенных дифференциальных уравнений / Д. Холл, Д. Уатт. - М.: Мир, 1979. 312 с.

13. Марчук, Г.И. Методы вычислительной математики. - М.: Наука, 1989. 608 с.

VELOCITY FIELD CALCULATION FOR VISCOELASTIC FLOW IN THE DEVICE WITH MIXER

© 2014 N.A. Gazizullin Kazan National Research Technological University

The laminar viscoelastic flow in the device with a mechanical mixer is carried out by the control volume method. The computational results are presented as streamline contours of secondary fluid circulation.

Keywords: device with mixer, mixing, viscoelastic fluid, circulation, streamlines

Nazym Gazizullin, Candidate of Technical Sciences, Associate Professor at the Higher Mathematics Department. E-mail: gnazym@gmail.com