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

МЕТОДЫ РАСЧЕТА ХИМИКО-ТЕХНОЛОГИЧЕСКИХ ПРОЦЕССОВ

УДК 532.517.2

Б. А. Снигерев, Ф. Х. Тазюков, М. А. Кутузова ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ ТЕЧЕНИЯ УПРУГОВЯЗКОЙ ЖИДКОСТИ В ВИСКОЗИМЕТРЕ С ПАДАЮЩИМ ГРУЗОМ

Рассмотрена модель течения полимерной жидкости в канале вискозиметра с падающим грузом. В качестве реологического конститутивного соотношения использованы: 1. модель упруговязкой жидкости Oldroyd-B; 2. ре-лаксационнаяя модель упруговязкой жидкости Фан-Тьен-Таннера (PTT), предсказывающая аномалию вязкости раствора полимера и наличие продольной вязкости в течениях, обусловленных нормальным напряжением. Для этих моделей получены численные решения задачи о падении груза при различных граничных условиях на стенках канала. Результаты расчетов удовлетворительно согласуются с численными и экспериментальными данными других авторов.

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

Согласно общей схеме течения в вискозиметре с падающим грузом, на груз сферической формы действуют следующие силы: сила Архимеда (Fbuoyancy); вес груза и суммарная сила вязкого трения (Fvisc). В случае движения груза сферической формы в вязкой жидкости, находящейся в неограниченном контейнере, суммарная сила вязкого трения Fvisc определяется в соответствии с формулой Стокса: Fvisc = 6nnrv. При достижении грузом состояния равномерного и прямолинейного движения v = const эти силы уравновешиваются.

Пусть груз сферической формы имеет радиус r и плотность р и движется под действием сил гравитации, трения и силы Архимеда в цилиндрическом контейнере радиуса R . Упруговязкая жидкость, находящаяся в контейнере, характеризуется плотностью р',

вязкостью при нулевой скорости сдвига По и характеристическим временем релаксации

напряжения. Полимерная жидкость может также содержать растворитель с вязкостью. п3 и обладать эффектом аномалии вязкости.

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

^иоуапоу +1^130 + (-тд) = 4яг3рд + 6пп™ - 4пг3рд = 0.

Откуда скорость падения груза V определяется по формуле

2 г2д .

V = -—(Р-Р).

9 п

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

Несмотря на кажущуюся простоту, течение жидкости вблизи груза достаточно сложное. А именно, течение вблизи передней и задней кромки груза происходит в основном под действием нормальных напряжений, а течение в канале между грузом и стенкой контейнера происходит под действием значительных сдвиговых напряжений. Относительная важность вклада различных факторов, влияющих на скорость падения груза, зависит в первую очередь от реологических свойств исследуемой жидкости, веса груза и относительных размеров груза и контейнера.

В настоящей работе приведено решение двумерной задачи об обтекании бесконечного цилиндра в плоском канале. Исследованы две схемы обтекания цилиндра, которые приведены на рис. 1.

Постановка краевой задачи

Рассмотрим постановку задачи течения упруговязкой жидкости в канале при обтекании цилиндра под действием перепада давления. Уравнения движения в терминах экстра напряжений имеют вид

д V — —

р (—+ V -V V) =-УР + V- т ; (1)

д I

V -V = 0 , (2)

где р - плотность жидкости, V - вектор скорости, Р - давление, т -девиатор напряжения.

Реологическое уравнение для данной модели имеет вид

ЛЧ^ Л**, Л**,

т = Т1 + т2 > (3)

где т1 - неньютоновская составляющая напряжения, определяемая по формуле

т + Х ( т1+ Ц тгб)) + Г2(Х)т1 = 2Л1 (О ). (4)

Ньютоновская составляющая общего напряжения т2 определяется из соотношения

Х2 = 2п 26. (5)

V

Здесь т1 - верхняя конвективная производная от тензора напряжения т1, п1 - вязкость

полимера, п2 - - вязкость растворителя, X - время релаксации напряжений.

Выбор параметров ^ , f2 определяет вид конститутивного реологического соотношения. При ^ = 0 и f2 = 0 получается контравариантная модель упруговязкой жидкости ОШгоуё-В, а выбор функций ^ _ в виде

~ ~ ~ <г< т)

^(т, О) = ^(Р- т + т-О), Г2 = е "" . (6)

определяет экспоненциальную модель упруговязкой жидкости Фан-Тьен-Таннера (РТТ). В уравнении (6) параметры £ртт и ^ являются материальными параметрами, определяющими поведение данной модели. Верхняя конвективная производная определяется как

V д т — — —

т = —+v-Vт-V V- т-т- (Vv)T , (7)

д I

а тензор скоростей деформации определяется соотношением

О = 1 (уу + ^т ). (8)

Уравнения (4)-(5) запишем в безразмерном виде введением следующих переменных

* х * у * и * V . Ut * ЬР * Ьт /Г1Ч

х = —, у =—, и = —, V = —, t =—, р =-------------, т =-----, (9)

1*1 и и 1_ п и п и

где и - характеристическая скорость задачи; 1_ - характерный размер, соответствующий ширине канала. Применение характерных масштабов (9) к уравнениям (1) - (5) приводит к безразмерному виду

V - -т

We ( т1+ МЗД) + Т2 (X) т = (1- в) (VV + VV ), (10)

Яе^+V-V V) =-vp +в V2 V, (11)

V-V = 0, (12)

где безразмерные параметры Яе и We являются числами Рейнольдса и Вайссенберга

риь х 1 и

Яе = риЬ We = —^, (13)

П Ь

а параметр в определяется как

в=-^. (14)

П1 +П2

Метод решения

Для решений уравнений (10)-(12) применяется метод конечных элементов. Для этого запишем «слабую» формулировку этих уравнений в виде [2,3,6]

— + V -V v)wdn = (‘(^Р + вv 2_ "~1

д t

| Яе( —^ + V -V v)wdn = | (^Р + вV2 V +Vт1)wdn ; (15)

|(V v)q dП = 0 ; (16)

п

|^е ( т1 + f1( т, О)) + f2 (X) т}w dП = | (1- в) (VV + VvT )w dП . (17)

п П

Далее область покрывается семейством четырехугольных элементов. Для скоростей и напряжений применяется квадратичная аппроксимация и линейная для давления. Интегралы на четырехугольнике вычисляются с помощью 4-х точечной квадратурной формулы Г аусса. После применения стандартной процедуры метода конечных элементов получаем матричные уравнения следующего вида

Re M vn+1 + в A vn+1 - D pn+1 = Re M vn - Re Cvn + D т п , (18)

D T vn+1 = 0 , (19)

где M - матрица масс; D - дискретизация оператора дивергенции; D T - транспонированная матрица оператора дивергенции; A - сеточный аналог оператора Лапласа; C - дискретизация конвективного члена.

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

WeM тп+1 = WeM тп - WeG(vn) тп + Fтп -С*(уп)т п, (20)

где M матрица масс; G -дискретизация оператора верхней конвективной производной; F - порожден дополнительными членами уравнения (17). Для стабилизации счета применяется дополнительная вязкость вдоль линии тока в виде

N * = N +а V-V N (21)

где для конвективных членов вводится весовая функция N , определяемая по формуле (21). Формулы для вычисления параметра а, зависящего от размера и средней скорости в данном элементе приведены в работе [2]. Для решения уравнений на каждом временном слое применяется прямой метод разложения Холецкого [7]. Решение начально-краевой задачи получено с помощью процедуры установления по времени. Для обеспечения сходимости решения по сетке задача решалась на последовательности сгущающихся сеток с числом узлов конечно-элементной сетки N = 12240, 24600.

Течение в канале, при обтекании цилиндра

Рассмотрим стационарное неизотермическое течение разбавленного полимера в канале, имеющего препятствие в виде цилиндра, схема которого представлена на рис. 1. Рассмотрены два варианта течения: - в первом, рассматривается обтекание цилиндра в канале, стенки которого имеют скорость и (рис. 1а) и второй, когда происходит пуазейлев-ское обтекание цилиндра в канале (рис. 1Ь). На рисунке 2 И и 2R обозначают ширину канала и диаметр цилиндра, х, у - декартова система координат.

Схема 1а

Рис. 1 - Схемы расчетных областей при обтекании цилиндра

В центре расчетной области х = у = 0 расположен цилиндр, входная граница расположена на расстоянии от центра Ц =15И, а выходная на 1_2 =30И . В силу этого, на входе задается полностью развитое течение полимера, определяющегося в безразмерном виде как

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

На твердых стенках канала ставятся условия прилипания жидкости и= V = 0, а на оси ставится условие симметрии V = 0, тху =0 . Для расчетов использовалась четырехугольная

конечно-элементная сетка, состоящая из N = 24600 узлов и М = 24100 элементов. Сетка конечных элементов существенно сгущается в области вблизи цилиндра для более точного учета больших градиентов напряжений.

Моделирование течения проводилось при следующих значениях параметров:

Re = 0.01, р= 0.11, £ = 0.3, £ = 1.0, We = 0.01, 0.5, 1.0, 3.0.

Реологическая конститутивная модель Фан-Тьен-Таннера

В этом разделе представлены результаты моделирования течения жидкости, удовлетворяющей конститутивному соотношению Фан-Тьен-Таннера (РТТ).

Одним из наиболее интересных и малоисследованных явлений реальных течений полимеров в канале при обтекании цилиндра является образование так называемого «отрицательного следа», возникающего только при движении груза в упруговязкой жидкости.

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

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

Результаты наших расчетов значений продольной составляющей скорости для схем течения, представленных на рис. 1, приведены на рис. 2.

Как следует из рис. 2, для схемы течения с подвижными стенками (рис. 2а) наблюдается максимальное значения скорости при We =0.5 и составляет и = 1.25, что означает 25% превышение скорости над значением скорости на оси канала. При дальнейшем увеличении времени релаксации наблюдается некоторое уменьшение пика скорости, но процесс замедления до скорости на оси канала происходит на большем расстоянии.

(22)

Эх д х д х д х д х

ди д Т д Т ху д Т ху д Т уу

^ ь хх =___________________ху. =___________ху. =___уу. = о р = о

Результаты численного моделирования

Аналогичное ситуация наблюдается в поведении профиля скорости в следе за цилиндром при обтекании в канале потоком Пуазейля (рис. 2Ь). Из рис. 2 также видно, что в случае близком к ньютоновскому поведению при малых значениях времени релаксации (We = 0.01) ситуация аналогична для обеих схем течения. Скорость монотонно восстанавливается из нулевого значения на цилиндре до скорости в основном потоке.

В дальнейшем определена такая важная характеристика течения упруговязкой жидкости, как разность главных напряжений а1 - а 2 =^№ + 4тХу , где Ы1 = т^ - туу. В соответствие с оптическим законом а1 - а2 = М(п1 - П2), где (п1 - П2) определяет главную разность коэффициентов преломления [9].

Рис. 2 - Профили скорости в следе за цилиндром, для модели жидкости РТТ при различных числах We

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

На рис. 3 показано распределение величины разности главных напряжений в следе за цилиндром для обеих схем обтекания цилиндра, указанных на рис. 1.

Из рис. 3 также следует, что при малых значениях числа Вайссенберга величина

о1 - о2 резко возрастает, а потом на расстоянии, равном 2 ширинам канала от задней кромки цилиндра, уменьшается до нуля. С ростом значений числа Вайссенберга, пик вели-

чины о1 - о2 снижается, однако релаксация главных напряжений происходит на гораздо большем расстоянии от задней кромки цилиндра (свыше 4 ширин канала).

В дальнейшем аналогичные исследования были проведены также для реологического конститутивного соотношения Oldroyd-B.

Реологическая конститутивная модель ОШгоу<^-В

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

Рис. 3 - Профили разности главных напряжений в следе за цилиндром

На рис. 4 показано распределение продольной компоненты скорости в следе цилиндра. Для схемы течения Пуазейля в канале (схема 2б) численные расчеты, в соответствие с рис.8Ь, показывают монотонный характер восстановления скорости. В случае течения с подвижными стенками (схема 2а) данный эффект тоже имеет место, только выражен он очень незначительно при We = 0.5. При увеличении числа We его проявление также исчезает. Таким образом, для обеих схем течения, представленных на рис. 1, эффект «отрицательного следа» незначителен или практически отсутствует. Данная ситуация свидетельствует о том, что на возникновение эффекта «отрицательного следа» большее влияние оказывает наличие свойства аномалии вязкости упруговязкой жидкости.

Далее проанализируем распределение разности главных напряжений. Как было показано, что в случае использования реологического соотношения РТТ, характеризующегося учетом аномалии вязкости, с ростом времени релаксации напряжений (ростом числа We), происходит некоторое уменьшение пика разности главных напряжений (рис. 3) с одновременным ростом расстояния на котором о1 - о2 релаксирует.

В случае использования реологического конститутивного соотношения Oldroyd-B ситуация меняется на противоположную. Как видно из рис. 5, с ростом времени релаксации напряжений в потоке (с ростом величины числа Вайссенберга) величина пика разно-

сти главных напряжений а1 - а 2 + 4тХу возрастает, а не убывает, как предсказыва-

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

Рис. 4 - Профили скорости в следе за цилиндром

Рис. 5 - Профили разности главных напряжений

Картины течения вблизи цилиндра

В дальнейшем рассмотрим картины течения, полученные для различных значений числа Вайссенберга. На рис. 6 обтекание цилиндра происходит слева направо. На рис. 6a, показаны линии уровня горизонтальной компоненты скорости и, на рис. 6б показаны изобары, на рис. 6c показаны изолинии сдвиговых напряжений т ху, на рис^ показаны линии уровня разности главных напряжений. Отметим, что пользуя^ возникновением эффекта

двойного лучепреломления, вызванного течением полимерной жидкости вблизи цилиндра, можно экспериментально наблюдать изолинии (П1 - П2), которые согласно оптическому закону пропорциональны разности главных напряжений о1 - о2 [9].

Зона максимальной величины разности главных напряжений (рис. бё) расположена вдоль поверхности цилиндра, причем наибольшие их значения расположены на боковых поверхностях, где наблюдаются наибольшие градиенты скорости. На кривых распределения скорости потока и отчетливо проявляется зона превышения скорости (рис.ба ). Также видно, увеличение концентрации напряжения на осевой линии, а вследствие этого и о1 - о2 в зоне за цилиндром. Отметим, что максимальные значения безразмерных напряжений при этом уменьшаются и их распределение по потоку становится более равномерным.

Рис. 6 - Картина течения жидкости РТТ для схемы 2а при We = 1.0

82

Как указывалось выше, реологическая модель упруговязкой жидкости ОШгоуё-В применяется для полимерных жидкостей (растворов или расплавов), обладающих постоянной сдвиговой вязкостью. Поэтому в случае использования данной модели проявление эффекта «отрицательного следа» за цилиндром практически не наблюдается. Небольшое увеличение скорости было обнаружено только для схемы течения (2а) при We = 0.5, однако при увеличении числа We данный эффект полностью исчезает.

Рис. ба, где показаны линии уровня скорости U, иллюстрирует данное явление. Также отметим, что расчеты по данной модели показывают увеличение значений напряжений вблизи цилиндра. Трудность расчета по данной модели состоит в том, что данная модель характеризуется, в частности, очень большими градиентами напряжений в тонком пограничном слое вдоль поверхности цилиндра. Значения напряжений в следе за цилиндром также монотонно повышаются при увеличении числа We . Это может быть объяснено тем, что данная модель не предсказывает свойство аномалии вязкости, и поэтому значение напряжений монотонно увеличиваются при увеличении числа We. Отметим, что расчеты по данной модели показывают также очень большие градиенты напряжений, как на поверхности цилиндра, так и на стенках канала. При этом значения напряжений гораздо больше, чем при расчете по модели РТТ. Поэтому расчеты полимеров по модели Олдрой-да являются хорошим тестом для апробации численных алгоримтов. Численные схемы пригодные для расчета нестационарных течений моделей данного типа будут гарантированно пригодны для моделей других типов вязкоупругих жидкостей.

Вычисление силы вязкого трения Ру|зС

Суммарная сила вязкого трения в проекции на продольную ось определяется по формуле

Fvisc = 0 [(-з +2П2 |Х+т„)СО$0 + (п {|Х+|у}) +Ч^е. (23)

Однако большинство исследователей вычисляют коэффициент коррекции Са, определяемый по формуле

Сс1 ' ^окБ = ^ или Са| = i .

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

Таблица 1 - Зависимость Са от значения числа We

We Са ОШгоу^В схема (2а) Са РТТ схема (2а) Са ОШгоу^В схема (2Ь) Са РТТ схема (2Ь)

0.0 102.12 102.76 132.36 132.36

0.5 89.46 79.25 126.41 120.96

1.0 81.96 70.45 124.12 118.96

3.0 - 66.45 - 110.78

Из таблицы 1 следует, что для заданных геометрических соотношений, интегральное влияние неньютоновских свойств приводит к уменьшению величины коэффициента коррекции. В результате расчетов коэффициента коррекции при обтекании цилиндров установлено, что свойство аномалии вязкости играет более значительную роль в снижении сопротивления цилиндра. Уменьшение сопротивления цилиндра по модели РТТ, учитывающей данное свойство является более значительным. Коэффициент сопротивления для плоского обтекания, когда поток ограничен двумя неподвижными стенками (схема 2б) превышает величину сопротивления в случае относительного движения стенок канала (схема течения 2a) приблизительно на 20%.

Заключение

Результаты численного моделирования стационарного обтекания цилиндра в канале с неподвижными стенками и в канале, стенки которого движутся с некоторой скоростью показали, что в обоих случаях вязкоупругая жидкость в следе за цилиндром проявляет нелинейное поведение, отличное от стоксового обтекания ньютоновской жидкостью. Данное свойство проявляется в образовании так называемого «отрицательного следа», профиль скорости за цилиндром носит немонотонный характер. Проведенное численное исследование данного эффекта показывает рост зоны «отрицательного следа2 при увеличении времени релаксации жидкости, что позволяет сделать предположение, что данный эффект связан с проявлением эффекта нормальных напряжений, а также аномалией вязкости вязкоупругого материала в сдвиговом течении в такой форме. Расчеты для модели жидкости Олдройд-Б, обнаруживают очень слабое проявление данного эффекта, что подтверждает утверждение о значительном влиянии в этом обтекании свойства аномалии вязкости вязкоупругого материала.

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

Литература

1. Хаппель Дж., Бреннер Г. Гидродинаимка при малых числах Рейнольдса. M.: Мир, 1976, 623c.

2. M.Aboubacar, J.P. Aguayo, P.M. Phillips, T.N. Phillips, H.R. Tamaddon, B.A. Snigerev, M.F. Webster. Modelling pom-pom type models with high-order finite volume schemes// J. Non-Newtonian Fluid Mech. 2005. V.126. P.207-220.

3. Ф.А. Гарифуллин, Б.А.Снигерев, Ф.Х.Тазюков, М.А.Кутузова, Амер Аль-Раваш. Численное моделирование обтекания цилиндра потоком упруговязкой жидкости Олдройда-Б //Тр. VII Между-нар. конф. “Нефтехимия -2005”. Нижнекамск. 2005. C 219-221.

4. Dou H.S., Phan-Tien N. Negative wake in the uniform flow past a cylinder // Rheol. Acta. 2003. V. 42. P.383-409.

5. Aboubacar M., Phillips T.N., Jahrom H.R., Snigerev B.A. High order finite volume methods for viscoelastic flow problems.// J.Comput. Physics. 1999. 16-40. 2004.

6. Полежаев В.И., Бунэ А.В., Верезуб Н.А. и др. Математическое моделирование конвективного тепломассообмена на основе уравнений Навье-Стокса. М.: Наука, 1987. 272 с.

7. Писсанецки С. Технология разреженных матриц. М.: Мир, 1988. 412 с.

8. Arigo M.T., McKinley G.H. An experimental investigation of negative wake behind spheres settling in a shear-thinning viscoelastic fluid.//Reol. Acta. 1998. V.37. P.307-327.

9. Виноградов Г.В., Малкин А.Я. Реология полимеров. М.: Химия, 1977. 440 с.

© Б. А. Снигерев - канд. техн. наук, ст. научн. сотр. ИММ КазНЦ РАН; Ф. Х. Тазюков - д-р. техн. наук, проф. каф. теоретической механики и сопротивления материалов КГТУ; М. А. Кутузова - асс. каф. процессов и аппаратов химической технологии НХТИ.