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

УДК 532.137.3

МОДЕЛИРОВАНИЕ ДВИЖЕНИЯ ВЯЗКОПЛАСТИЧНОЙ СРЕДЫ В ЦИЛИНДРЕ, СОВЕРШАЮЩЕМ ЗАТУХАЮЩИЕ КРУТИЛЬНЫЕ КОЛЕБАНИЯ

ГЛ. Вяткин, И.В. Елюхина

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

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

Математическая формулировка задачи. Кратко напомним основные условия движения [1]. Пусть пустой цилиндр подвешен вдоль своей оси на упругой нити и совершает вокруг нее крутильные колебания с периодом г0 и декрементом затухания ¿>0 . При заполнении его вязко-пластичной средой изосинхронный режим колебаний в общем случае нарушается, и колебательный процесс можно охарактеризовать параметрами т- 2АТт и 3 - 21п|ах /а21, где АТт - разница между двумя соседними моментами времени, когда угловое смещение цилиндра а обращается в нуль; а},а2 - соседние экстремальные значения а (¡«я^ > |#2|). Период г>г0 вследствие того, что при увлечении жидкости стенками цилиндра возрастает эффективный момент инерции подвесной системы; д > <50 в связи с дополнительной диссипацией механической энергии, обусловленной вязким трением в жидкости.

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

- уравнение движения жидкости

дТ £ 9

- уравнение движения цилиндра

^ + ^ + а = (2)

йт я ат

- начально-краевые условия для (1), (2):

А г/

= 0; (3)

£/(£0) = 0, = и(09Т) = 0, а(0) = а0~6\ —

аТ дТ

7=0

Вяткин Г.П., Елюхина И.В. Моделирование движения вязкопластичной среды _в цилиндре, совершающем затухающие крутильные колебания

- реологическое уравнение состояния вязкопластичной среды (модель Бингама - см. [2, 3, 4] и пр.)

I1 + Вг//Е>) при ° -

О^ = 0 при сг < Вт\

(4)

где

ди и А МЯ2 л 50 . % = —- — = Вт- 0

$ » Н ^ 2К урск

и - , Т = £0 =11/<1, £ = г/с1, й =у/у/д0 , =2я/т0\ (5)

Л - отношение момента инерции жидкости в вискозиметре (МЛ2 / 2) и пустой подвесной системы (К ); Вт - число Бингама; й - толщина пограничного слоя; О - безразмерный второй инвариант тензора скоростей деформации, для рассматриваемых условий течения В =| |;

- безразмерная ая компонента тензора скоростей деформации; М - масса среды; Р -момент сил, приложенных к цилиндру со стороны среды; - циклическая частота колебаний пустого цилиндра; К - внутренний радиус цилиндра; г - радиальная координата (г = О на оси цилиндра); / - время; V - азимутальная компонента скорости; а0 - начальное угловое смещение цилиндра; <р - угловая координата; V , р, 50 - кинематическая вязкость, плотность и предел текучести среды соответственно; а - безразмерный второй инвариант тензора напряжений; а^ - безразмерная £<р~я компонента тензора напряжений; задача (1)-(4) решается численно.

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

При заполнении вискозиметра вязкопластичной средой около оси цилиндра всегда присутствует твердое ядро, где сдвиговые напряжения не превосходят предела текучести. В потоке также имеется тонкий твердотельный слой, возникающий у поверхности цилиндра, перемещающийся затем к ядру, граница которого движется в это время от центра, и при достижении ядра сливающийся с ним; в следующую четверть периода происходит обратное движение; далее в новом полупериоде при £ = вновь образуется такой слой (в общем случае их может быть несколько) и т.д. Зависимость, к примеру, радиуса ядра от времени является периодической функцией; колебания £ = обычно нарастающие, а кривая, отражающая среднее значение, является вогнутой возрастающей функцией. Эти колебания не гармонические, даже если не учитывать резкий скачок при слиянии ядра и слоя, а спектр = в целом аналогичен таковому для

а = а(Т) (основная частота = равна Фг/г1, где г1 - соответствует гармонике с наибольшей интенсивностью на спектре а = а(Т)). В застойных, или твердотельных, зонах скорость II по координате £ изменяется линейно: й\] / - V / ¿, = 0, и, в частности, в случае, когда твердое ядро заполняет весь объем вискозиметра, и{^,Т)-да/с/Т • £. Таким зонам соответствуют прямолинейные участки, начиная от £ = 0, на профилях скорости сдвига и скорости среды, а также участки с искривлением профиля скорости сдвига при смене знака И.

С ростом номера колебания N в процессе затухания область «твердотельного» течения растет, и при N > Nтв твердое ядро заполняет весь объем вискозиметра. В этом случае эффективный момент инерции системы достигает своей наибольшего значения + К), а вместе с ним наибольшим становится и период, прямо пропорциональный корню квадратному из этой величины. Значение 6 тогда минимально и совпадает с <50 ввиду отсутствия диссипации меха-

нической энергии вследствие вязкого трения. Зависимость безразмерной частоты колебаний заполненного жидкостью цилиндра Л = т/т0 от N является монотонно возрастающей и тем быстрее достигает значения Лтв (при N > N тв режим колебаний так же изосинхронный, как и для ньютоновской среды), чем больше Вт . Зависимость д - д{Щ сначала растет до некоторого значения 8т, одинакового для всех чисел Вт при прочих равных условиях, а затем сравнительно быстро падает до д ~ д0. При увеличении числа Вт кривая 5 = смещается к оси ординат и становится монотонно убывающей. В подобной ситуации наблюдаемые в эксперименте эффекты, связанные с вязкопластическим поведением среды, обусловлены, прежде всего, не гармоническим законом колебаний даже в отсутствии переходных процессов. Поэтому для большей наблюдаемости таких эффектов должна быть большей разница между значениями частоты Л1 и декремента дх в начале колебательного процесса и значениями Лтв и дт,

соответственно. К тому же, желательно, чтобы колебания затухали дольше, чем достигались значения д0 и Лтв, число N было достаточным, но в то же время чувствительности дЛ/дВт и дб/дВт (а также и дЛ/ Э/У, дS/дN) имели высокие значения. Выбор оптимальных условий эксперимента можно выполнить аналогично приведенному, например, в [5].

Способы численного моделирования. Расчет по идеальной модели для вязкопластичной среды, так же как и учет конечной длины цилиндра, значительно усложняет численные формулировки и в общем случае не позволяет обеспечить требуемую точность и в ряде случаев сходимость к истинному устойчивому решению при заданных требованиях к эксперименту. Так, в частности, численные модели, основанные на определении на каждом временном слое радиуса твердого ядра вязкопластического течения аналогично выполняемому при стационарных течениях, к примеру, в капиллярной реометрии, не позволяют корректно промоделировать твердотельные прослойки, наличие которых изменяет напряжение на стенке цилиндра, внося зачастую существенные ошибки в закон колебаний. К тому же, при определенных условиях эксперимента и само определение радиуса ядра вызывает значительные трудности ввиду неоднократной смены знака в распределении скорости сдвига вдоль радиуса цилиндра. Заметим также, что наличие слагаемых в системе уравнений цилиндра и среды, связанных с торцевыми эффектами (необходимых в случае не 'длинного' цилиндра), делает овраг функции качества, являющейся критерием соответствия расчетных и экспериментальных данных, на плоскости (у,а0) еще более пологим и практически невозможным корректную оценку свойств среды при реализуемой на практике точности измеряемых параметров (подробнее о построении функции качества, планировании оптимального эксперимента и пр. - см., к примеру, [5]).

В связи с вышесказанным для описания идеального бингамовского поведения в качестве основной рабочей модели была выбрана модель Ы-лчвсоБку (см, например, [6. 7, 8]), обеспечивающая при решении этой внутренней гидродинамической задачи хорошую сходимость с идеальной вязкопластичной жидкостью. Напомним, что в модели движение нетекучей компоненты рассматривается как движение ньютоновской среды, вязкость которой намного превосходит вязкость текучей компоненты:

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

ка ~103... ка ~104 расхождение в рассчитанных значениях параметров колебаний составляло даже в наихудшем случае менее 0,1 %.

иуи и с- ,

при В < В0;

при В > О,

(6)

альному вязкопластичному поведению значение ка принималось достаточно высоким ~103 (обычно ка~ 102). Такое приближение дает удовлетворительные результаты: при

Вяткин ГЛ., Елюхина И.В.

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

II. Неидеальная вязкопластичная среда. Модели вязкопластичного поведения. При моделировании не идеального бингамовского поведения рассмотрим все возможные варианты:

1) модель Ы^БСОБку (6) с различными значениями модельного коэффициента ка;

2) экспоненциальную зависимость эффективной вязкости от модуля скорости сдвига, предложенную в [9], и также, как и Ы^бсоб^у, широко используемую на практике (см., например, [10, 11] и пр.)

= |1 + ^[1-ехр(-тО)]|^; (7)

где т - параметр, контролирующий экспоненциальный рост напряжения, который обычно в расчетах при моделировании бингамовского поведения принимается ~102;

3) реологическое уравнение состояния нелинейно вязкой среды с показателем степенного реологического закона п < 1 (модель Оствальда - Вейля - см., например, [12])

а4<р=Ьйп-10^ (8)

где Ъ - / , к - постоянная степенного реологического закона.

Для описания нелинейностей, возникающих на экспериментальных кривых течения бинга-мовских сред в диапазоне малых скоростей сдвига, используются модели нелинейного вязкопластичного поведения: Кессона, Балкли-Гершеля и пр., а также модель Шульмана, сочетающая в разных нелинейных мерах пластичность и нелинейную вязкость и обобщающую большинство известных моделей. Особенности колебательных процессов при заполнении вискозиметра такими средами (на примере реологической модели Балкли-Гершеля) были изучены в [1].

Далее примем следующие условия эксперимента: = 10, А = 0,1, 50 < 0,001, и будем моделировать течение среды с числом Бингама Вт = 0.5. Для наглядности реологических схем расчета приведем также кривую течения - зависимость напряжения от модуля скорости сдвига (см. рис. 1), показывающую как приближает та или иная модель поведение реальной среды к идеальному вязкопластичному в диапазоне О, реализуемом в вискозиметре. Параметры моделей (6)-(8), принятые для расчетов, следующие:

- модель (6) (пунктирные линии на рис. 1 и далее): 1 - кг = 10, 2 - кг = 100 ;

- модель (7) (штрихпунктирные линии): 1 - т -10, 2 - т = 100;

- модель (8) (сплошные линии): 1 - Ь = 0,75 и п = 0,9 , 2 - Ъ = 1,25 и п — 0,7 .

Рис. 1. Кривые течения неидеальных вязкопластичных сред

Изменение параметров колебаний при движении цилиндра. Зависимость параметров колебаний от их номера для сред со сложными реологическими свойствами продемонстрирована на рис. 2. В случае идеального вязкопластичного материала после установления 'твердотельного5

режима течения (когда N > 7Уте) значения 6 ~ 30 и Л ~ Хтв = (1 + А)112 (повторим, что результаты расчетов по идеальной среде получены в рамках модели Ы^зсозку с коэффициентом А^-103

и приведены в [1]). В остальных расчетных ситуациях при N > Nтв еще присутствует диссипация вследствие вязкого трения, и 8 > 80; а Л < Лтв, т.к. не вся жидкость внутри цилиндра выполняет роль «присоединенной» массы. Из рис. 2 можно увидеть какой степенью неньютоново-сти обладают среды с рассматриваемыми реологическими уравнениями состояния согласно этой особенности. Значения параметров колебаний при аналогичных условиях эксперимента для ньютоновской среды А,ньют -1,014 и §ньют -0,069, а Хтв -1,049. Наилучшее приближение к идеальной бингамовской среде из представленных на рис. 1 вариантов обеспечивает модель Ы-у1зсо8ку с ка =100, что согласуется результатами, касающимися изменения параметров колебаний во времени (рис. 2а, б).

Рис. 2. Зависимость параметров колебаний от их номера

На зависимостях 8 ~8{N) и Л~Л{Ы) (рис. 2а, б) имеются участки, где 8, Л -const. Они, например, в модели bi-viscosity, отвечают вискозиметру, полностью заполненному нетекучей средой. В терминах модели (6) при низких значениях ка это участки, которым соответствует течение ньютоновской среды с вязкостью, равной ка. Так, при вязкости, равной ка =10 (вместо 1) параметры колебаний при расчете по зависимости для ньютоновской среды X -1,041 и 5 - 0,0856 (соответствуют кривой течения (прямой), обозначенной знаком *, на рис. 1), При вязкости, равной 100, значения X -1,0486 и 5-0,0113, что отвечает величинам, указанным на рис. 2 а, б. Выход на режим с S(N), Л(Л0 - const в случае проведенного расчета происходит быстрее при более низком ка. Это можно объяснить тем обстоятельством, что при меньшем ка интервал значений D, при которых вязкость постоянна, шире и к тому же в этом случае декремент затухания выше, что обеспечивает более быстрый выход при меньшем N . Участки д, Л - const

горизонтальны (к примеру, при ка =10), т.е. как и в случае колебаний идеальной среды при N > Nme, но значения 80 и Лте в этом случае не достигаются, т.к. такой вариант является слабым приближением идеальной бингамовской жидкости. Смещение кривых д = S(N) и Л = Л(Ы) с изменением модельного коэффициента в некотором смысле аналогично реализуемому при изменении числа Вт и прочих равных условиях для идеальных сред, т.е. чем сильнее влияние пластичности, тем больше смещается к оси ординат, а затем и вовсе исчезает, максимум функции 8 -8{N); для сред с большими модельными коэффициентами значения Л1 и дЛ/дВт выше и пр. Так, к примеру, с ростом ка (при ка~102 и выше) происходит подобное смещение указанных кривых, и выход на режим с постоянными значениями частоты и декремента затухания обычно происходит быстрее.

При увеличении чисел Бингама при некотором Вт = Вт' напряжение в среде становится недостаточным для сдвига, обеспечивающего возникновение текучей фазы, и с самого начала коле-

Вяткин ГЛ., Елюхина И.В.

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

бательных процессов вискозиметр полностью заполнен твердой зоной. При таком сдвиге, хотя и малом даже для ТУ = 1 (с ростом N при одном и том же режиме течения амплитуда В = В(Т) становится меньше), ввиду перестройки напряжений в потоке возможно существование течения и при N > 1. При малых Вт с ростом этого числа кривые 8 = 8(Ы) и А = А^) смещаются последовательно, а при Вт происходит скачок и далее при всех Вт параметры колебаний для идеальной вязкопластичной среды 8 ~ 80 и А~ Лтв, При уменьшении Вт значения Л1 и 81 приближаются к Лньют и 8ньют, а при Вт<~ 0,01 они практически им равны (по рис. 2 можно проследить, как влияет наличие пластичных свойств на 8Х). Зависимости параметров колебаний от их номера в этом случае также аналогичны приведенному на рис. 2 а, б: декремент сначала растет, а потом падает; период растет монотонно, начиная от Лньют, и в процессе затухания колебаний реализуется больший интервал значений Л и 8 .

Для ньютоновских сред (рис. 3 в [13]) величина Л уменьшается с ростом в основном на интервале е (2... 12), а поведение 8 зависит от величины : при > с ростом декремент уменьшается, а при < - растет, где - 4,2. Эти закономерности подтверждаются и в случае вязкопластичных сред. Так, представим последние средами с очень высокой вязкостью (в ка раз выше, чем у ньютоновской - см. (6)) и в качестве рассмотрим величину

ы - I ^¡^' г#е к о- ~ усредненная по полупериоду кажущаяся вязкость. В процессе колебаний кп растет, и значение £оы уменьшается. Значение -4,2, и, к примеру, соответствует экстремуму 8-8{И) на рис. 2а. Заметим, что если с самого начала колебательных процессов реализуется значение < 4,2 , то зависимость 8 = 8(Ы) - монотонно убывающая.

Поясним эту качественную картину с иной стороны. Рассмотрим бингамовскую среду как псевдопластичную с реологическими свойствами, соответствующими (8). Тогда

£опу ~ , где В ~ усредненное по полупериоду значение В при £ = . С ростом N

амплитудное значение скорости сдвига уменьшается, и кажущаяся вязкость ЪВП~Х растет, а значение £0|1У падает. Для ньютоновских сред при значениях —> 0 с уменьшением период растет слабо, а декремент при 80 -0 устремляется к нулю [13]. Также и для нелинейно вязких (псевдопластичных) сред в результате колебаний происходит выход на асимптотический режим, характеризуемый мало изменяющимися во времени значениями Лас. При уменьшении в интервале е (3.. .1) происходит сильное падение 8 , что соответствует участку кривой 8 - 8(Ы) после максимума; при —> 0 значения декремента 8ас соответствуют параметрам колебаний сильно вязкой среды.

Некоторые особенности течения вязкопластичных сред в цилиндре. С помощью реологической модели нелинейно вязкой среды можно наглядно объяснить распределение скорости течения вдоль радиуса цилиндра. Для таких сред можно выделить два вида течения: при вязкости больше и меньше ньютоновской, в зависимости от чего граница области развитого течения приближается или удаляется от цилиндра по сравнению с ньютоновской. При п < 1 кривые течения имеют выпуклость вверх, к оси напряжений (рис. 1), при п> 1 - вниз, а для ньютоновской среды, при п = 1, это прямые. Вязкость (кажущаяся, определяемая отношением напряжения и скорости сдвига) при малых скоростях сдвига в случае п < 1 имеет значение выше, чем для среды с п = 1, до точки пересечения кривой, ей соответствующей, с прямой, отвечающей п = 1, а с ростом В

становится ниже. Глубина проникновения пропорциональна кажущейся вязкости ЬВп~~1, которая

для псевдопластичных сред при данных условиях эксперимента (когда В < 1) падает с ростом п и В. Это означает, что для дилагантных сред (с п > 1) область развитого течения находится вблизи стенки цилиндра и при численном моделировании расчеты необходимо проводить в интервале £е ' где определяется как £/(()...£,7*) ~ 0, а для псевдопластичных - в ряде

случаев брать большее число точек у оси цилиндра. Для бингамовской жидкости (как и для псевдопластичной) также характерен рост области развитого течения по сравнению с ньютоновской, что хорошо можно наблюдать, к примеру, при колебаниях тела в безграничной среде. Это обусловлено, прежде всего, увлечением жидкости твердыми прослойками. Заметим, что если напряжение сдвига между «твердотельным» слоем и ядром недостаточно для возникновения течения в этой области, то при N, близком к Nme, «твердотельная» фаза может быть представлена только ядром (с присоединенным к нему слоем).

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

Литература

1. Елюхина И.В., Вяткин Г.П., Бескачко В.П. Новые возможности крутильно-колеба-тельного метода Швидковского Е.Г.: идентификация реологической принадлежности среды// Вестник ЮУрГУ. Серия «Математика, физика, химия». - 2003. - Вып. 3. - № 6(22). - С. 108-115.

2. Bird R.B., Dai G.C., Yarusso В J. The rheology and flow of viscoplastic materials // Rev. in Chem. Eng. - 1983. -№ l.-P. 1-70.

3. Barnes H.A. The yield stress - a review or «panta roi» - everything flows?// J. Non-Newtonian Fluid Mech. - 1999. - № 81. - P. 133-178.

4. Shelukhin V.V. Bingham viscoplastic as a limit of non-newtonian fluids// J. Math. Fluid Mech. -2002.-№4.-P. 109-127.

5. Вяткин Г.П., Елюхина И.В. Разработка методов параметрической идентификации сред Оствальда-Вейля по результатам вибрационной реометрии// Вестник ЮУрГУ. Серия «Металлургия». - 2004. - Вып. 4. - № 8(37). - С. 22-27.

6. Keentok М., Milthorpe J.F., O'Donovan Е. On the shearing zone around rotating vanes in plastic liquids: theory and experiment // J. Non-Newtonian Fluid Mech. - 1985. - № 17. - P.23-35.

7. Alkhatib M.A.M, Wilson S.D.R The development of Poiseuille flow of a yield-stress fluid // J. Non-Newtonian Fluid. Mech. - 2001.-№ 100.-P. 1-8.

8. Beverly C.R., Tanner R.I. Numerical analysis of three-dimensional Bingham plastic flow // J. Non-Newtonian Fluid Mech. - 1992. - V. 42. - P. 85-115.

9. Papanastasiou T.C. Flows of materials with yield // J. Rheol. - 1987. - № 31. - P. 385-404.

10. SmyrnaiosD.N., Tsamopoulos J.A. Squeeze flow of Bingham plastics// J. Non-Newtonian Fluid Mech. - 2001. - № 100.-P. 165-190.

11. Blackery J., Mitsoulis E. Creeping motion of a sphere in tubes filled with a Bingham plastic material // J. Non-Newtonian Fluid Mech. - 1997. - № 70. - P. 59-77.

12. Fang P., Manglik R.M., Jog M.A. Characteristics of laminar viscous shear-thinning fluid flows in eccentric annular channels // J. Non-Newt. Fluid Mech. - 1999. - № 84. - P. 1-17.

13.Швидковский Е.Г. Некоторые вопросы вязкости расплавленных металлов. - М.: ГИТТЛ, 1955.-206 с.

Поступила в редакцию 26 января 2005 г.