УДК 532.137.3
ОСОБЕННОСТИ ТЕЧЕНИЯ НЕНЬЮТОНОВСКИХ СРЕД В ОСЦИЛЛИРУЮЩИХ РЕОМЕТРИЧЕСКИХ СИСТЕМАХ
И.В. Елюхина
Путем численного моделирования экспериментов изучены особенности течения нелинейно вязких и вязкопластичных сред в осциллирующих рео-метрических системах на примерах крутильно-колебательного вискозиметра в режиме затухающих колебаний и вибрационного вискозиметра в режиме вынужденных колебаний.
Введение
При планировании корректных экспериментов по идентификации реологической принадлежности жидкостей и определения их параметров необходимо учитывать особенности течения сред различных реологических типов. Рассмотрим подробнее два типичных нестационарных случая: затухающие колебания крутильно-колебательного вискозиметра (рис. 1) и вынужденные колебания вибрационного вискозиметра (рис. 2) (см., например, [1]).
В подобных периодических течениях на практике обычно реализуются случаи, когда полная деформация остается сколь угодно малой, хотя ее мгновенная скорость может быть высокой. В этом случае состояние простой жидкости с затухающей памятью можно охарактеризовать как линейное вязкоупругое, а течение описать единственной материальной функцией - комплексной вязкостью. В этом случае для установившегося режима вискозиметрические уравнения могут быть получены аналитически, и в детальном изучении локальных особенностей течений таких сред нет необходимости. Поэтому далее остановимся на рассмотрении течений неньютоновских жидкостей, в реологические уравнения состояния которых входят нелинейно вязкие или вязкопластические составляющие.
Математическая формулировка задачи
Крутильно-колебательный метод
Математическую модель эксперимента для случая длинного цилиндра в цилиндрической системе координат представим в следующем виде:
1) уравнение движения жидкости
ди = д^£<р 2^£оО дТ ~ д£ + £
2) уравнение движения цилиндра \\^\ 1
'/А
Рис. 1. Схема метода:
1 - упругая нить; 2 - цилиндр, совершающий затухающие крутильные колебания вокруг своей оси; 3 - исследуемая жидкость
Рис. 2.Схема метода:
1 - упругий элемент; 2 - пластина,
совершающая плоские колебания под действием гармонич. силы;
3 - исследуемая жидкость
(1)
d a dT‘
■ + 2 А 0
dT
■ + a = P.
3) начально-краевые условия для (1), (2):
da Иг/
U(£,0) = О, U(£T) = —£0, U(О,T) = 0, a(0) = a0~60, ^ dT dT
= 0.
4) реологическое уравнение состояния
4.1) для ньютоновской среды
<5р = ,
4.2) для нелинейно вязкой среды по модели Оствальда-Вейля (см., например, [2])
ар= ЬЯрЯт-\
4.3) для вязкопластической среды по модели Бингама (см., например, [3])
[(1 + Вт/^р при Я > Б0,
при Я < Я0,
kvDfr
(2)
(3)
(4)
(5)
(6)
где
4
P = —г Ag Л £0^ Л1#=£
dU U Л MR2
Dpm =---------------, A =------
Л d£ £ 2K
Bm = -
vpq0
b =
„m-і ту-q0 Kv
vp
D0 =
Bm
kG-l і
U = V/dq0 , T = q01, £0 = R / d , £ = r / d , d = ^jv / q0 , q0 = 2п/т0, A 0 = S0 /2n; (7)
t - время, « - угловое смещение цилиндра, ск0 - начальное смещение, R - внутренний радиус цилиндра, d - толщина пограничного слоя, Kv и га - постоянная и показатель степенного реологического закона, Вт - число Бингама, <70 - предел текучести, v - кинематическая вязкость среды, р - плотность среды, V(r, t) - азимутальная компонента скорости , т0 - период собственных установившихся колебаний, q0 - циклическая частота колебаний пустого цилиндра, S0 -логарифмический декремент затухания собственных установившихся колебаний, Dp - £р -я компонента D, D - второй инвариант D, D - тензор скоростей деформации, Ор - р -я компонента тензора напряжений; P - момент сил, приложенных к цилиндру со стороны среды, A -отношение моментов инерции среды в цилиндре 0,5MR2 и пустой подвесной системы K относительно оси цилиндра, M - масса среды, величины P, M и K отнесены к единице длины цилиндра. Для рассматриваемого течения D = |Dp|. При моделировании вязкопластического поведения в (6) принята модель bi-viscosity, в которой для лучшего соответствия модели истинному вязкопластичному поведению модельный коэффициент ко ~ 103. Затуханием колебаний в отсутствие среды далее пренебрегаем: 80 = 0.
Вибрационный метод
Математическую модель вискозиметрических экспериментов представим как
ди д2и
дТ
д z
dT2
+ х = sin—+ 2A—=-
і
dz
dx / dT
t =0
= 0 , x(0) = 0, U(z,0) = 0 , U(0,T) = y • dx/dT, U(<*>,T) = 0 ;
(8)
(9)
(10)
для среды Оствальда-Вейля:
kf = mb
dU / d z
ktr = b
dU / d:
z=0
(11)
T=0
k
tr ?
z=0
m-1
m-1
Расчет и конструирование
для среды Бингама:
kf = К = ka - при
dU / д z
< Dn
kf = 1, ktr = | 1 + Bm / dU / д z
(12)
где x = xk / F .
T = m01,
X = m0/ m.
m0 = k / m0,
при dU / dz > D,
z = z / d , d = -yJv/m0 , A = Svpm0 / F
b = mm-Kv /(vp), Bm = <r0 /(vpm0), у = F /(kd), U = V /(dm0), V = Vx - скорость колеблющейся пластины, m0 - масса подвесной системы, k - жесткость упругого элемента, m и F - частота и амплитуда вынуждающей силы, S - площадь поверхности пластины, x - линейное смещение пластины, z - координата вдоль оси Z, z = 0 - на пластине; в = 42Ay, где в - см. [2]; описание остальных параметров - см. после (7).
Сопряженные задачи о движении зонда и среды (1)-(3), (5) или (6) и (8)-(10), (11) или (12) решены численно. Для решения системы нелинейных уравнений в работе использован метод прямых. Полученная система обыкновенных дифференциальных уравнений интегрировалась, в частности, методом Рунге-Кутта четвертого порядка с контролем точности и автоматическим выбором шага по времени, методом Адамса пятого порядка точности в форме Нордсика и пр. Для интегрирования жестких систем, возникающих при определенных условиях эксперимента, был использован метод Гира шестого порядка точности. Производные по координате аппроксимировались разностными отношениями с пятью узловыми точками, обеспечивающими точность порядка четвертой степени шага по координате.
Результаты и обсуждение
Нелинейно вязкие среды
Рассмотрим простейший тип осциллирующего течения - на пластине. Для нелинейно вязких сред можно выделить два вида течения: при вязкости больше и меньше ньютоновской, в зависимости от чего граница области развитого течения приближается или удаляется от пластины по сравнению с ньютоновской (рис. 3 - для дилатантной среды). Это объясняется тем, что глубина проникновения пропорциональна кажущейся вязкости bDm~1, определяемой отношением напряжения и скорости сдвига, которая для дилатантных сред при данных условиях эксперимента (когда D < 1) падает с ростом m и при уменьшении D. Так, рис. 3б соответствует малым D (DqI < 0.1, где D'0 = D'|^ , D' - скорость сдвига). Кривая течения для m = 2 на рис. 3а проходит
ниже прямой для ньютоновская среды с m = 1, и граница области, где U ~ 0, находится ближе к пластине для среды с m = 2. При в = const с ростом у, т.е. с ростом, например, амплитуды вынуждающей силы F , расширяется интервал значений скорости сдвига, проходимых в процессе колебаний, и при |D'0| >> 1 граница находится существенно дальше от пластины, чем для ньютоновской жидкости. На практике это может означать, что предположение о безграничности среды уже не допустимо и влиянием стенок пренебрегать нельзя.
2
-т = 2 -т = 1 1/
//
/ у /
0 0.5
1 IA,
а)
б)
Рис. 3. К зависимости и =и(£) для дилатантной среды при малых скоростях сдвига
На частотных спектрах напряжения и скорости сдвига на пластине (рис. 4), построенных для режима установившихся колебаний, выявлено появление нечетных гармоник вынуждающей силы, обнаруженных ранее в [4] в упрощенной модели с сосредоточенными параметрами. Распределение интенсивностей пиков спектра, в т.ч. вид их огибающей, определяется нелинейными свойствами жидкости и условиями эксперимента, что позволяет исследовать свойства таких сред в рамках Фурье реологии.
а)
b = 1, в = 1, т = 3, 1 = 1,1, у = 0,01 б)
Рис. 4. Спектр (а) и зависимость скорости сдвига (б) на пластине от времени
Данные особенности наблюдаются и при течении в крутильно-колебательном вискозиметре. Заметим, что в т.ч. при малых скоростях сдвига, реализуемых в режиме затухающих колебаний, среды с m > 1 относятся к дилатантным, а с 0 < m < 1 - к псевдопластичным. При численном моделировании с учетом того, что для дилатантных сред область развитого течения находится вблизи стенки цилиндра, расчеты необходимо проводить в интервале £е [£,£], где £ определяется как U(0... £, T) ~ 0 , а для псевдопластичных - в ряде случаев брать большее число точек у оси цилиндра.
Вязкопластичные среды
При заполнении крутильно-колебательного вискозиметра вязкопластичной средой около оси цилиндра всегда присутствует твердое ядро (зона 1 на рис. 5), где сдвиговые напряжения не превосходят предел текучести. Также в потоке имеется тонкая твердотельная прослойка (зона 2 на рис. 5), возникающая у поверхности цилиндра, перемещающаяся в процессе колебания к ядру, граница которого движется в это время от центра, и при достижении ядра сливающаяся с ним (рис. 5-7). В таких зонах скорость U по координате £ изменяется линейно: dU/d£ - U/ £ = 0, и, в частности, в случае развитого по всему сечению твердотельного течения U(£,T) = da/dT ■ £ . Застойным зонам соответствуют прямолинейные участки, начиная от £ = 0, на рис. 6, 7 (зона 1) и искривление профиля скорости сдвига при смене знака D' на рис. 7 (зона
2); это области с |D'| < D0 на рис. 5.
Расчет по идеальной модели для вязкопластичной среды (т.е. когда Op = Dp = 0 при
0р| < о0 в (6) и пр.), так же как и учет конечной длины цилиндра, значительно усложняет численные формулировки, в т.ч., существенно повышая время расчета, и в то же время не обеспечивает требуемой точности и в ряде случаев сходимости к истинному решению. Так, например, численные модели, основанные на определении на каждом временном слое радиуса твердого ядра вязкопластического течения аналогично выполняемому при стационарных течениях, к примеру, в капиллярной реометрии, не позволяют корректно промоделировать твердотельные прослойки, наличие которых изменяет напряжение на стенке цилиндра, внося часто существенные ошибки в закон колебаний. Использование модели bi-viscosity здесь дает удовлетворительные результаты: при kO~103... kO~104расхождение в рассчитанных значениях составляет менее 0,1 %.
Расчет и конструирование
0.1-с/ D'
D
О
0.02
0.015
0.01
0.005
0
0.005
' V ■
*
■ ■'
/ л / /2
1 /
1
1
U
0.3
0.2
0.1
0
-0.1
-0.2
1 -
2 - 3
4 - ,. -
\
\ \ \ 1
1
D
0.067
0.033
0
0.033
0.067
-0.1
\
\
\ \ \ '' \
\ \ \ \
. \ \ \
10
11
: 9 10
Рис. 5. Зависимости
U=U(& и l)'=l(Ç)
- застойные зоны
Рис. 7. Зависимость D=D(£)
0 3 6 9
Рис. 6. Зависимость и=и(£)
£0 = 12, А = 0,2, Вт = 0,4;
Кривые 1-4 соответствуют моментам времени в течение 1-й четверти по-лупериода (крив. 1, 2) и 2-й четверти (крив. 3, 4)
При течении в вибрационном вискозиметре при значениях г, прилегающих к правой границе расчетной области, образуется неподвижная твердотельная зона, а вязкопластическое течение перемежается твердотельными прослойками (обычно 1-й - 2-мя). Для такой жидкости характерен рост области развитого течения по сравнению с ньютоновской (т.е. смещения расчетной границы и (тс,Т) = 0 в (10)) вследствие увлечения жидкости твердыми прослойками, что качественно согласуется с результатами для псевдопластичных сред.
Заключение
В работе обсуждены особенности, на которые следует обращать внимание при численном моделировании экспериментов по идентификации реологической принадлежности и свойств неньютоновских сред - нелинейно вязких (модель Оствальда-Вейля) и вязкопластичных (модель Бингама). Так,
1) выявлен рост области развитого течения для вязкопластичных и псевдопластичных жидкостей и уменьшение этой области для дилатантных сред;
2) проанализирована динамика развития твердотельных зон в процессе затухания крутильных колебаний,
3) установлено появление нечетных гармоник вынуждающей силы на спектрах напряжения и скорости сдвига на пластине;
4) выполнены рекомендации по численным расчетам таких течений в осциллирующих рео-метрических системах.
Литература
1. Елюхина И.В., Вяткин Г.П. Параметрическая идентификация нелинейно вязких свойств жидкостей вибрационным методом затухающих колебаний//Вестник ЮУрГУ. Серия «Машиностроение». - 2005. - Вып. 6. - № 1 (41). - С. 6-11.
2. Balmforth N.J., Craster R.V. A consistent thin-layer theory for Bingham plastics// J. Non-Newtonian Fluid Mech., 1999. - № 84. - Р. 65-81.
3. Fang P., Manglik R.M., Jog M.A. Characteristics of laminar viscous shear-thinning fluid flows in eccentric annular channels// J. Non-Newtonian Fluid Mech., 1999. - № 84. - Р. 1-17.
4. Wilhelm M., Maring D., Spiess H.-W. Fourier-transform rheology//Rheol. Acta, 1998. - № 37. -Р. 399-405.