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

ТЕХНОЛОГИЯ ПОЛУЧЕНИЯ И УТИЛИЗАЦИИ ЭНЕРГЕТИЧЕСКИХ КОМПОНЕНТОВ

УДК 532.519.6

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

A.A. ВАХРУШЕВ, A.M. ЛИПАНОВ, A.B. ВАХРУШЕВ Институт прикладной механики УрО РАН, Ижевск, Россия

АННОТАЦИЯ. Рассмотрена задача моделирования двумерных течений неньютоновских жидкостей в криволинейных каналах. Проведены численные расчеты параметров псевдопластической среды, описываемой с помощью реологической модели Пауэлла - Эйринга, и исследованы физические явления, возникающие при локальном изменении температуры стенки канала.

ВВЕДЕНИЕ

В процессе движения жидких сред в различных устройствах, таких как трубопроводы, перерабатывающее оборудование, при извлечении вязких ископаемых из несущих пластов зачастую геометрия каналов, по которым происходит движение, является определяющим фактором, влияющим на основные параметры потока. Наличие выступов, изгибов приводит к возникновению вихревых течений, уменьшению расхода жидкости вследствие возрастающего градиента давления. Дополнительными факторами, усложняющими анализ указанных физических процессов, являются реологические характеристики жидкости: вид реологического закона, определяющий, каким образом реагирует среда с той или иной внутренней структурой на внешнюю нагрузку; пределы изменения коэффициента вязкости, его связь с температурным полем и т.д. Многофакторное исследование комплексного влияния указанных выше факторов на параметры течения обусловливает актуальность использование математического моделирования. Настоящее исследование является развитием работы [1], в которой рассматривалось течение неньютоновских жидкостей в прямом канале при наличии структурных (дискретных) изменений вязких свойств среды под действием температурного поля. Рассмотрена задача течения в плоском двумерном канале, состоящем из прямых участков, соединенных изогнутым ("коленом"), в предположении, что на входе в изогнутую часть на одной из адиабатических стенок происходит нарушение теплоизоляции, и жидкость локально охлаждается.

ПОСТАНОВКА ЗАДАЧИ

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

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

Общая форма уравнений Навье -Стокса с учетом переменной вязкости может быть представлена в виде [7]:

х

_ 1 ър 1

(1)

р • г д(р г д(р г2 дг

Т»г

\ др ар^) 1 с>(г

п

(2)

р дг г д(р г дг г

1 дг ■ V --^ +-:

г д(р дг

= 0,

(3)

где Я ,5ГГ,8 ,5 - компоненты тензора напряжений.

+

г дер дг г

8п. = 2?] ^ г/ I о

<3г

(4)

Если предположить, что вязкость жидкости ц не меняется по времени и внутри области интегрирования, то уравнения (1),(2) упрощаются до вида:

5у

ч>

дч.

1 дР

V ЭУ У V — + уг—2. + -*-—=---

дI дг г дер г р- г дер

+ П

\ 2 дуг

Ду._ - -4- +

<р

2 дер

(5)

где

дvr дуг i V V 1 ар —Л + —- + —------- =---

д( дг г дер г р дг

Ау 2

2

Г

г2 дер

А =

15 1

+--+ —г-

дг" г дг г

9 '

(6)

(7)

Переход к безразмерным переменным при решении задачи течения жидкости в криволинейном канале позволяет записать (5)-(6)

£>у

<р _

дР 1 +

Ж дер Яе

г~ г~ дер

(8)

Л

дР 1 +

дг Яе

г 2 2 а

г г дер

(9)

где = ^

число Рейнольдса, характеризующее тип течения (ламинар-

ное/турбулентное); й,, /?2 - внутренний и внешний радиусы кривизны канала; 1/0 - характерная скорость потока, втекающего в канал (в частности, берется максимальное значение скорости на оси канала, соответствующее профилю Пуазейля); 77, р - кинематическая вязкость жидкости и плотность среды, соответственно.

Уравнение неразрывности (3) для использования в дальнейших преобразованиях можно представить как:

Л

—- + V г + г —- = О дер дг

с учетом условия г ф 0.

Уравнения энергии движения (8)-(10) записаны для физических переменных "скорость-давление". В первом приближении будем считать, что расход жидкости на входе в канал постоянный. Тогда из расчетов можно исключить давление и свести исходную систему уравнений к двум неизвестным - "вихрь-функция тока".

Для вывода уравнения переноса "вихря" со использовались соотношения:

со

поляр.

1

г

д гу

<р

дг

£4

дер

(П)

V

д\р

поляр.

<р

дг

уг =— г

1 ду/

поляр.

дер

(12)

со

поляр.

5У поляр. (

дг2 г дг

поляр

- +

1 д2у/

поляр. _

дер'

= АI//

поляр.

(13)

Легко показать, что переменные бУполяр, ,<//(|ошф в пределе, при радиусе кри-

визны г —» оо, совпадают со своими аналогами в декартовых координатах г/, V, со, \[/ . Это утверждение помогает осуществить переход на границах прямого и изогнутого участков трубы от полярных к декартовым координатам. Далее опустим нижний индекс и будем обозначать "вихрь" и "функцию тока" как со и , соответственно.

Для получения уравнения переноса вихря в полярных координатах исключим давление из уравнений закона сохранения импульса (5)-(6).э используя представление со в виде (11):

V дсо дсо V

—--+ уг — = —

г дер дг г

+у.

1 <4 1 а2

V.

+

г дер дердг г дер'

+

V

ч>

+ -

1 д* д\п 1 эу„ 1 52у.

+

<р

о О

г" г дг дг" г" дер г дердг

_ <р

у Зу у <Э2у у д2 <р <р + _!?_ ч> V и

у.

г дер г дердг г дер~

у у

<р г

ду

+

<р

+ у

д2\ д\>.. у. 32у.

+

г г дг дг~ г дер г дердг

(И)

Для проведений численных расчетов важной является недивергентная форма уравнения переноса. Она позволяет использовать метод характеристик и численные схемы, использованные автором при расчете параметров течения в прямом канале [1].

Было получено уравнение переноса "вихря" со в недивергентной форме:

Ра>_д<о 1 д(уа>) | д(у,-а>)

£)/ д{ г дер дг

(15)

Полученное уравнение (15) не учитывает диссипативные силы, связанные, в случае течения вязко -пластической жидкости, с силами вязкого трения. Для получения источниковых членов в нем были использованы уравнения (1)-(2), в предположении, что вязкость жидкости // не меняется по пространству и времени. Это приводит (для безразмерных координат) к формуле

Ф

Яе

д'со 1 дсо 1 д'со

-Г +--+ ~-7

дг г дг г' д<р~

= ДСУ.

(16)

Выражения для тензора напряжений (4) в цилиндрических координатах использовались для получения диссипативной части уравнения переноса "вихря" в виде:

Ф(0 = Д(т/ - со) + ф„.

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

Таким образом, получена диссипативная составляющая уравнения переноса в цилиндрических координатах:

0> =

(О

и N 4 а>

МП—т-

г~ о (рог

1 д?] д27]

\

г дер дердг ^

2 (Уу

г дг2

/

\

>

дг/ + 1 д"т] дг г дер2

2

г дг1

1 д2у/ ду/ г дер" дг

(17)

/

Необходимо отметить, что знаки перед слагаемыми в формуле (17) могут меняться, в зависимости от предполагаемого положительного направления "завихренности" со.

Для расчета параметров течения в изогнутом канале, с учетом реологической модели Пауэлла - Эйринга, необходимо записать уравнение энергии в полярных координатах. Воспользуемся общей формулировкой уравнения энергии для полярных координат [7]:

дт V дТ дт \дд \д( ,л ~ 1

<Э/ г дер дг г дер г дг

ду

<р

[дер

+ V

+

д_ дг

v

р

ч г у

1 0У. + -

г дер

дг

(18)

1 1 дТ

1 дТ

- компоненты вектора плотности энергии (в со-

Ре г дер Ре дг

ответствии с законом Фурье).

С учетом соотношений для тензора напряжений (4) и переходя к недивергентной форме записи уравнения энергии, получим

д( г дер дг г Ре

(19)

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

Ф7. = 2

ду

<р

дер

+ v.

+

'дО2

ч дг у

+

а

дг

г.. \

v

<р

v г у

+

г дер

(20)

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

РЕЗУЛЬТАТЫ РАСЧЕТОВ

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

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

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

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

Конечно-разностный аналог уравнения переноса "вихря" со запишется следующим образом:

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

где

1

ао ' % + аи • œM / "аи' ' œi-i j - af '60

(22)

Ц.со

J_ и a r

(23)

a r;A<p

j[a>MJ-2 соц+й)^

(24)

Цсо

_ - iy--

/'¡Аг

'••_, 1 г 1

(25)

В представленных зависимостях (22) и (23) величины a{j и Д; - характеристики

компонент скорости ( v^ и vr, соответственно) через грани контрольного объема [5,6]. С

учетом данных характеристик, частные производные по пространству, стоящие в конвективной части уравнения переноса, вычисляются с учетом направления движения жидкости. Это делает смысл конечно-разностного уравнения (21) более физическим, по сравнению с использованием центральных разностей при вычислении частных произ-

Bv.F dvrF

водных и -fcr .

Для решения системы линейных алгебраических уравнений, получаемой при использовании численной схемы (21), воспользуемся процедурой ADI (Alternative Directions Implicit) [2]. В соответствии с данным методом, разобьем задач}' на три временных слоя: n,n + j,n + 1,-атакже разделим тангенциальное и радиальное направления:

»4 At _,

"--¿"¿У

о Re *

>1+

о =С°'»

(26)

ûÇ] + At L\(Û

At T-) w+i

--L:œ

« Re У

= co 2.

У

(27)

В результате получим две трехдиагональные системы уравнений вида

(28)

а] ■ <oiJ+l + Ъ) ■ Cùtl + с] ■ а> х =fj,

где

(29)

а? =

Ъ? =

СГ =

/Г =

А( |2 1

а: -

М

г;Аср '' Яе г2Агр

■>»

. Д/ ,, ,2ч 2 1 +-(«,;'-«:-) +

М

г^ср

А1 21

ал -

А/

(30)

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

Ь-

А/ |2 _ 1 А/

Ру

1 А/

Дг

= 1 +

Ке А г2 2 • Яе гЛг '

Лср(Р" Р'-> Яе г2А(р2 '

(31)

Д/ „21 1 А/

-—В1.-Аг

+

1 А/

Яе Аг 2-ЯегЛг'

Л =

со0. \

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

Полученные трехдиагональные системы конечно-разностных уравнений (28),(29) решаются методом "прогонки" [5,6], описанным ранее.

Уравнение Пуассона (13), связывающее "завихренность" со с "функцией тока" в полярных координатах, можно представить в конечно-разностном виде, используя центральные разности. Данные аппроксимации частных производных дают второй порядок точности по пространству.

1

г.Аг

1

1

.^ ~]+ ^Т- 2' ^у + П-1 ]+ "2 • ^ + Ум.,] = ■(32)

Решение системы линейных алгебраических уравнений выполнено методом последовательной верхней релаксации [2,5,6]. Для этого, предварительно сгруппировав соответствующие слагаемые, приведем равенство (32) к виду

у*

1

2 • (/?;+!)

1 г,'к(р

где коэффициенты а п /?у вычисляются как а - = — гу. • Дг/>, р/ =--.

Сложим правую часть равенства (33) с тождеством уц - ул = 0 и перегруппируем слагаемые

1

улу = Уи + 2 (/?2 + 1)' + + + «,) • + (/?; - а,) • - 4а>, -2(^+1)- щ

+

(34)

В итоге, предполагая, что в ячейках с индексами / -1 и / -1 значения известны на новом п +1 временном слое, получим

у Г = К +

2-(Л2+1)

+ (/?,2 - а,) • - 4а*< - +1) • ^

+

(35)

где т - релаксационный параметр, для которого выполняется условие: 0 < т < 2.

Неявная схема обладает лучшей сходимостью при решении уравнений, описывающих движение сред со сложной реологией (рис.3).

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

Вязкость Шш&Ш

Рис 3. Реологическая поверхность, описывающая поведение псевдопластической жидкости

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

Для расчета течения по составной трубе потребовалась стыковка уравнения переноса "вихря" на границе участков. Сравним формулы для "завихренности" со в декартовых и полярных координатах:

_ди__ду_ дуг

^декарт. ~0у ~дх ' ^ " а7 + 7 Г ' дер '

По размерности обе величины совпадают, а с учетом того, что указанные части канала стыкуются под прямым углом, на граничной области у^ = и, у;. = у . Также, при выборе достаточно маленьких пространственных шагов расчетной сетки Ах, Ау, АериАг в пределе будут выполняться соотношения г • А ер —> Ах, Ау = А г . Тогда очевидно, что основной вклад в выражении Аа> = ¿Уполяр - <Удскарт., основываясь на (36), будет составлять у(р/г . Введем обозначение новой функции со

^ = где = -у-, (37)

которой соответствует новое уравнение переноса

дсо 1 5у со ду-со 1 д~?1'СО 1 дп-ео д"п-со — +---5-+ —-= —---— +---'—- + +—'-т-

9/ г дер дг г~ д г д дг~

где диссипативная функция Ф^ вычисляется с помощью соотношения (17), а С- определяется как

дб 1 ду -8 ду-8 — +----— +——

д( г дер дг

1 д2?]-3 дт]-8 [ д2?]-5)

кг2 д г д дг2

(39)

С помощью формул (37)-(39) проведена серия расчетов, как с учетом постоянной вязкости среды (рис.5), так и при наличии зоны охлаждения и повышении вязких свойств жидкости (рис.6).

На основе экспериментальных данных известно [8], что при течении вязкой жидкости по составному каналу (рис.2), возникают вихревые течения. Их форма, направление движения и области образования отображены на рис.4.

Рис. 4. Возникновение вихревых течений при движении жидкости по изогнутому каналу

8.3185 0.1224 -0.0657

а

Рис. 5. Течение в изогнутом канале при постоянной вязкости

0.6780 0.5092 0.3404 0.1715 0.0027 -0.1661 -0.3350 -0.5038

Рис. 6. Течение в изогнутом канале при наличии скачка вязкости

Данное явление наблюдается и в численном эксперименте (рис.5): обратное течение возникает на входе в изогнутую часть трубы и на некотором расстоянии от выхода жидкости из данного участка. Области обведены на рисунке.

Для расчета течения вязкой жидкости в изогнутом канале границы расчетной области принимались достаточно удаленными (порядка 5-6 поперечных размеров канала) от изогнутой части. Это условие использовалось для того, чтобы вихревые зоны, образующиеся в потоке, не оказывали существенного влияния на выходные границы канала.

Расчеты проводились так же, как и в задаче течения жидкости по прямому каналу, при числе Рейнольдса Яе = \0 в невозмущенном потоке, а в локальных областях вязкость 71 увеличивалась до 50 раз.

В случае постоянной вязкости (рис.5) поток жидкости, входящий по прямому участку канала в изогнутый, ''прилипает" к нижней стенке. В результате, в верхней час-

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

Наибольшие энергетические потери, связанные с вязким трением, наблюдаются в потоке в областях, отмеченных сплошными линиями на графике распределения "завихренности" со (рис.5(с)).

Течения с локальными областями изменения вязкости в изогнутом канале обладают еще более сложной структурой (см. рис. 6). На графике (Ь) данного рисунка показано распределение вязкости в случае охлаждения участка нижней твердой границы на входе в "колено". Встречая вязкую "преграду", поток отклоняется к верней стенке. Также как и в случае с течением в прямом канале, наблюдается явление "биения струи". Эффект данного явления усиливается появлением третьей вихревой зоны (на выходе из изогнутого участка). Она представляет собой развившуюся циркуляционную область при течении жидкости с постоянной вязкостью за счет дополнительных "колебательных" движений основного ядра потока.

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

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

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

A.A. ВАХРУШЕВ, A.M. ЛИПАНОВ, A.B. ВАХРУШЕВ

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

1. Липанов A.M., Вахрушев A.A., Вахрушев A.B. Математическое моделирование течения жидкости с переменной структурой // Химическая физика и мезоскопия. -2005. Т7. N1. - С.23-30.

2. Андерсон Д., Таннехилл Дж., Плетчер Р. Вычислительная гидромеханика и теплообмен: В 2т. Т.2 - М.: Мир, 1990. - 392с.

3. Флетчер К. Вычислительные методы в динамике жидкостей. В 2-х т. - М.: Мир, 1991. -Т.1 - 502с. Т.2-552с.

4. Лойцянский Л.Г. Механика жидкости и газа. - М.: Наука, 1973. - 848с.

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

6. Роуч П. Вычислительная гидродинамика. - М.: Мир, 1980. - 616с.

7. Слеттери Дж. С. Теория переноса импульса, энергии и массы в сплошных средах. -М.: Энергия, 1978.-448с.

8. Повх И.Л. Техническая гидродинамика. 2-е изд. - Л.: Машиностроение, 1976. -504 с.

SUMMARY. Simulation of non-Newtonian 2-D flows in bended channels has been considered. Numerical calculations of pseudo-plastic medium parameters were carried out using Paule - Eyring rheological model. Physical processes into account occurring due to local change of channel wall temperature were taken.