CC BY
23
УДК 536.2
Ф. Х. Тазюков, Ф. А. Гарифуллин, Э. Р. Кутузова
ЗАДАЧА О ТЕЧЕНИИ ВЯЗКОУПРУГОЙ ЖИДКОСТИ ЧЕРЕЗ ПЛОСКОЕ 8 : 1 СУЖЕНИЕ
Ключевые слова: метод контрольных объемов, модель ЕЕМЕ-Р.
Целью настоящей статьи является описание постановки задачи изотермического ламинарного сходящегося течения жидкости в канале с резким симметричным сужением 8:1. В качестве реологического конститутивного соотношения была выбрана модель ЕЕМЕ-Р. Численное решение проводилось методом контрольных объемов с использованием программной среды ОрепЕоат.
Keywords: control volume method, FENE-P model.
The purpose of the present paper is a description of the isothermal laminar converging liquid flow problem formulation through the axisymmetric planar 8:1 contraction. The FENE-P model was chosen as a rheologial constitutive relation. Numerical method was conducted using a control volume method and an OpenFoam software environment.
Введение
При производстве полимеров особый интерес представляет изучение поведения течения вяз-коупругих жидкостей, а именно входные потоки и потоки в устройствах с особой геометрией типа сужение, расширение, сужение-расширение и т.д. Численное изучение задач плоского течения имеет большое значение в силу выявления новых свойств течений, ранее смоделированных для определенного набора параметров. Трудность изучения их заключается в сложности моделирования и предсказания течения жидкостей при конечных значениях времени релаксации. Данная проблема представляет особый интерес для ученых по всему миру, начиная с 1987 года, когда эталонной проблемой была поставлена задача плоского течения с коэффициентом сужения 4:1. Для упруго-вязких жидкостей наблюдается эффект возникновения циркуляционного течения вблизи острой кромки канала, именуемый в литературе как lip vortex, и его объединение с угловым циркуляционным течением. В настоящей статье рассматривается плоское осесимметричное стационарное течение с коэффициентом сужения 8:1. Детально изучается поведение lip vortex, его взаимодействие с угловым течением при различных числах We.
Реологическое соотношение FENE-P
В качестве механической модели молекулы принимается гибкая гантель, представляющая две бусинки, соединенные между собой пружинкой. Движение гибкой гантели в потоке растворителя запишется в следующем виде:
2(e) -.(d) -,(Ь)
-(e)
2_
где Р; - сила упругости пружинки, соединяющая бусинки;
р( )=-Ь|(?ги0-1 Г|) - сила трения Стокса, действующая на бусинки гибкой гантели, возни-
кающая вследствие обтекания молекулы вязким растворителем. Сила трения Стокса (сила лобового сопротивления) действует на сферические объекты с очень маленькими числами Рейнольдса.
L¡ - коэффициент трения для i-той бусинки, определяемый как L¡=6nnR
r¡ - скорость /-той бусинки
Но - скорость среды, обтекающей бусинку
L - транспонированный тензор от градиента скорости
L=(Vu)T
r¡ - радиус-вектор -той бусинки;
F¡ - случайная сила, возникающая в результате теплового (броуновского) движения молекул растворителя и их столкновений с бусинками гибкой гантели.
Пренебрегая инерционными эффектами, содержащими массы бусинок, уравнение движения запишется в следующем виде:
А —> ^ 1 -(e) -(Ь)
r¡=Uo+(Lr¡)+-(F¡ +F¡ )
L¡
В дальнейшем пользуемся следующими гипотезами:
1) Полимерная жидкость состоит из вязкого растворителя и растворенного в нем полимера
2) Раствор является разбавленным. Это означает, что взаимодействие макромолекул полимера в растворе модно пренебречь.
3) Полимерные молекулы представляются в виде гибких гантелей.
4) Ориентация макромолекулы и ее растяжение полностью описывается вектором конфигурации Q =r2 —r-t, соединяющим концы гибкой гантели.
5) В масштабе длины гантели течение является однородным.
6) Вследствие многочисленных столкновений с молекулами растворителя, вектор конфигурации также является случайным вектором и подчиняется функции распределения L(Q, t).
Броуновские силы являются случайными величинами, то в соответствии с определением энтропии S=-k ln L и свободной энергии FCB=-TS, где
T - абсолютная температура, k - постоянная Больц-мана.
—»(Ь) FCB д F: = —CB =-k^-TinL
dr=
дп
(Ь)
Сила имеет чисто энтропийную природу, а работа этой силы на перемещение ^ равна
d г и расходуется на изменение свободной энер-(Ь)
гии dFсв, то есть dFсв=F¡ dr.
С учетом этих формул дифференциальное уравнение для каждой бусинки в отдельности: dr1 _ ^ кТ д ?1 dt 1 _ дг— _
dr2 кТ д ?2
0= —2-- г2+--—!ш
dt 2
1 д^2
Вычитая второе уравнение из первого, получим:
dr1 dr2 _ „ - „ кТ д кТ д ?2
dt dt 21 _ дг1 _ дг2 _2 _
или
dc:i___ д
—=-С-кТ_12 — !п_-_12? dt 12 дСС 12
11
где ^2= — + —, 1Г=1Г1=-1::2
1 2
Подставляя преобразованное дифференциальное уравнение движение бусинок в уравнение непрерывности функции распределения конфигурации макромолекул Ь(СД), получим:
dL
"дГ
(LQ)
дь
д2
-(e)N
=kTLi2 —+ Li2—(L F )
до д^ до
В этом уравнении второе слагаемое определяет влияние гидродинамических взаимодействий. Третье слагаемое показывает влияние сил, возникающие под воздействием броуновского движения. Четвертое слагаемое демонстрирует взаимодействие между бусинками в модели, образующееся за счет силы упругости пружинки.
Умножая уравнение на величину СС и проведя процедуру осреднения по ансамблю, получим следующее уравнение:
^ ----__ ~Т —>—>(е)
^ <СС> -ЦСО - <СО> - =2кТ_12 <01? >
или
2(e),
(СС) =2кТЬ121-2Ь12 <<31? >
V
где () - верхняя конвективная производная, удовлетворяющая принципу материальной объективности, оператор (..) означает процедуру осреднения.
В дальнейшем воспользуемся выражением Крамерса для связи компонент тензора напряжений
с компонентами тензора (С? ):
(е)
Т=-пкТ !+пзу+п <С? >,
где первое слагаемое является шаровым тензором и определяет вклад гидростатического давления в общее напряжение, второе - вклад вязкого трения в растворителе и связан с ньютоновским реологическим конститутивным соотношением, третье - оп-
ределяет вклад упругости пружинок в суммарное напряжение.
Проанализируем третье слагаемое. Рассмотрим куб, имеющий грань, ориентация которой определяется направлением вектора нормали V. Известно, что число пружинок, имеющих вектор кон-
формации в интервале от С! до С+dС , равно —> 3
п-1_(О^С , где распределение вероятности ЦС) уже не является равновесным распределением вероятности.
Вектор конформации С! выделенной пружинки, проходящей через выделенную грань куба, определяется ненулевым скалярным произведением вектора нормали V на вектор С!, то есть V-С!. Таким образом, величина (V-о) • (п-_(О^С3) является усредненным числом пружинок, имеющих размеры вектора конформации в интервале от С! до Q+dQ и пересекающих воображаемую выделенную грань.
Сила, возникающая в результате суммарного действия гипотетических пружинок величины (С!, СМ С!) на грани поверхности куба, равна
2(e)
(v q) • (n-L(Q)dQ3) •
(e)
где сила ^ соответствует силе упругости изолированной пружинки.
Суммируя по всем возможным длинам вектора конформации, получим результирующую силу упругости
-, Г -> —>(е) —> 3 _ —> —>(е)
?=Уп- I (0-1= )_(С^С3^п <0-1= >
С другой стороны, суммарная сила ? связана с напряжением Тр, действующим на единичной грани выделенного куба, ориентированной в соответствии с направлением вектора нормали. Эта связь определяется соотношением
?=^Тр
Исходя из последних двух уравнений, полимерная част тензора напряжений определяется по
~ -> —>(е)
формуле Тр=п {СО Р1 >.
Уравнение Крамерса может быть переписано в форме Гизекуса
п V
12
сзсг
-.-.(e)
Выделяя величину n {CD F1 > из уравнения V
(e)
<QQ> =2kTLi2l-2Li2 <QF >
получим:
_ —»(e) ni-;^)
n <QF > =--Si+nkTI
2 12
Подставляя полученную формула в выражение Крамерса, получим:
2
Конкретный вид компонент тензора напряжений зависит только от задания конкретного вида силы упругости.
Закон взаимодействия представлен в виде закона Уорнера (рис. 1).
• продольная вязкость - в продольных течениях зависимость продольной вязкости от скорости продольного деформирования [1].
Математическая постановка задачи
Изотермические течения неньютоновских жидкостей в сужающихся каналах описываются уравнениями движения и неразрывности[2]:
Этот закон обладает следующими свойствами:
а) при малых деформациях пружинки зависимость между силой и растяжением остается практически линейной, то есть близкой к закону Гука ;
в) дальнейшее увеличение деформации пружинки приводит к увеличению ее жесткости;
с) пружинка не может растягиваться более, чем на величину .
где р - плотность жидкости, - вектор скорости, - давление, - девиатор напряжения.
В соответствии с принципом расщепления напряжений девиатор напряжения представляется как совокупность неньютоновской и ньютоновской составляющих:
Здесь - вязкость полимерной составляющей жидкости, - вязкость растворителя, - характерное время релаксации, - безразмерный параметр, характеризующий степень растяжения данной макромолекулы, - тензор конфигурации, - тензор скоростей деформации.
Уравнения движения в декартовой системе координат:
Рис. 1 - Зависимость силы взаимодействия между бусинками от деформации
Реологическое конститутивное соотношение запишется в следующем виде:
Уравнения эволюции тензора конфигурации в проекциях на оси координат:
где Л- время релаксации, - кинетическая энергия броуновского движения молекул, - число гибких гантелей в единице объема или числовая плотность гантелей, - кинетическая энергия.
Реологическое соотношение БЕКЕ-Р предсказывает следующие свойства жидкости:
• аномалия вязкости - зависимость вязкости от сдвиговой скорости;
• упругие свойства - конечное время релаксации напряжений;
dv du
A
'xy
Axx+Ayy
3L2
Уравнения, связывающие компоненты полимерной части тензора напряжений с компонентами тензора конфигурации:
'/ л , \
т = П1
Txx л
т = П1
Tyy л
Ax
1
V
A +A
xx yy
1
1
3L2
У
A
yy
V1
Axx+Ayy
1-4
T = П1
Txy л
3L2 A
xy
A +A
xx yy
V1
3L2
В результате приведения исходных уравнений к безразмерному виду можно получить, что в уравнения движения входят числа Рейнольдса, Вайссенберга и коэффициент ретардации [3] ре= ЕН-, We=р=-^
п -
П1+П2
Рис. 2 - Схема канала с резким сужением 8:1
Граничные условия.
На входе в канал задаются скорости:
v=0, u=const, p=const На выходе из канала задается условие установившегося течения:
Для обеспечения условия установившегося течения длина узкой части канала принималась равной 15Hi.
-=0, v=0, p=0
dx ' ^
На твердых стенках канала задается условие прилипания:
v=0,u=0
На оси задается условие симметрии: du
— =0,v=0 dy 1
Метод контрольных объемов
Одним из наиболее часто используемых в инженерных приложениях является метод контрольных объемов. Метод обладает свойством консервативности и достаточно просто в интерпретации. Расчетную область разбивают на некоторое число непересекающихся контрольных объемов таким образом, что каждая узловая точка содержится в одном контрольном объеме. Дифференциальное уравнение интегрируют по каждому контрольному объему. Для вычисления интегралов используют кусочные профили, которые описывают изменение искомой функции между узловыми точками. В результате находят дискретный аналог дифференциального уравнения, в который входят значения искомой функции в нескольких узловых точках.
Полученный подобным образом дискретный аналог выражает закон сохранения искомой функции для конечного контрольного объема точно так же, как дифференциальное уравнение выражает закон сохранения для бесконечно малого контрольного объема. Одним из важных свойств метода контрольного объема является то, что в нем заложено точное интегральное сохранение таких величин, как масса, количество движения и энергия на любой группе контрольных объемов и, следовательно, на всей расчетной области. Это свойство проявляется при любом числе узловых точек, а не только в предельном случае очень большого их числа. Таким образом, даже решение на грубой сетке удовлетворяет точным интегральным балансам. [4]
Расчеты проводились с помощью метода контрольных объемов на неравномерной сетке 30x90 со сгущением вблизи острой кромки канала 300:1(рис. 3).
Рис. 3 - Неравномерная сетка расчетной области
В тестовом режиме были проделаны расчеты и на других сетках: более густых и более разреженных. Однако, результаты расчетов для более густых сеток показали те же результаты, что и явилось обоснованием использованием выбранной сетки. Расчеты проводились для различных чисел
Вайссенберга ^е= в интервале от 0 до 900 при значении Рейнольдса (Ре= —) Ре=0.01 и коэффи-
циента ретардации
(ß= ) ß=1,
V ls+Пр/ 9
где ns
вязкость
растворителя, пр- вязкость полимера.
Задача изначально решалась как нестационарная, стационарное решение получилось методом установления при ^
L
1
Выводы
Литература
В настоящей статье представлена математическая постановка задачи о течении вязкоупругой жидкости модели БЕКЕ-Р в канале с плоским сужением 8:1. Дано описание, вывод реологического конститутивного соотношения БЕКЕ-Р с указанием гипотез, лежащих в основе модели, и свойств. Для решения поставленной задачи использовалась программная среда ОрепБоат и метод контрольных объемов.
1. Кутузов А.Г..Дисс. докт техн. наук, Казанский Государственный Технологический Университет, Казань, 2010, 351 с.
2. Ф. А. Гарифуллин, Ф. Х. Тазюков, М. А. Ахмадиев, Р. С. Шайхетдинова, Вестник Казанского технологического университета, 17, 7 14-17(2014)
3. Ф.Х. Тазюков, Х.А. Халаф, К.М. Алиев, Р.С. Шайхет-динова, Вестник Казанского технологического университета, 15, 4, 113-115(2012)
4. S. Patankar, Numerical heat transfer and fluid flow, Hemisphere Publishing Corporation, New York, 1980, C. 7475, 80-82.
© Ф. Х. Тазюков - д.т.н., профессор каф. ТМиСМ, КНИТУ, tazyukov@mail.ru; Ф. А. Гарифуллин - д.т.н., профессор каф. ТКМ, КНИТУ, rami21@yandex.ru; Э. Р. Кутузова - ассистент каф. АССОИ, КНИТУ, elvira.kutuzova@list.ru.
© F. Kh. Tazuykov - Ph.D, Professor of the Theoretical Mechanics and Strength Materials department,KNRTU, tazyukov@mail.ru; F. A. Garifullin - Ph.D, Professor of the Techology of Engeneering Materials department, KNRTU, rami21@yandex.ru; E. R. Kutuzova - assistant of the Automated Data Acquisition and Processing Systems department, elvira.kutuzova@list.ru
