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

Вычислительные технологии

Том 14, № 2, 2009

Моделирование турбулентного перемешивания однородной жидкости методом крупных вихрей

У. С. АБДИБЕКОВ, Д. Б. ЖАКЕБАЕВ, Б. Т. ЖУМАГУЛОВ Казахский национальный университет им. аль-Фараби, Алматы e-mail: [email protected], [email protected]

Приводятся результаты моделирования гидродинамики перемешивания однородной жидкости в цилиндрической области, внутри которой расположен пропеллер. Лопасти пропеллера вращаются под углами а = 90°, а = 45° к основанию цилиндра. Построена математическая модель, составлен численный алгоритм, получено решение задачи. Использованы осредненные по ансамблю уравнения Навье— Стокса в цилиндрических координатах.

Ключевые слова: гидродинамика больших чисел Рейнольдса, численное моделирование, турбулентное перемешивание, метод крупных вихрей.

Введение

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

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

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

2

пастей пропеллера, расположенного на высоте -Н. Рассмотрены случаи, когда лопасти

3

© ИВТ СО РАН, 2009.

находятся под углами 90 и 45° к основанию цилиндра. Предлагаются математическое описание данной проблемы и реализация численного решения задачи.

1. Основные уравнения

Основными уравнениями гидродинамики для описания турбулентного перемешивания являются осредненные уравнения Навье—Стокса, записанные в цилиндрической системе координат (г, е, г):

(дьт дьт

Р{-Ж + ^ -дГ

ье дьт

+--— + ^ я

дг

1 д

+ г дг V дг

дьт

+

г де 1 д2уг

+

дьт г д 2ьт

г2 де2 дг2

др дг

2 дье

г2 де 1

1 д 1 д д

(гТтт ) + -ТГ" (Тте) + ТГ" (Тт х )--(Ттт )

1 г дг г де дг г

2

V

е

/дье дье ье дье дье ьтье\ 1 др Р \ дЬ т дг г де х дг г ) г де

1 —(г—\ + + д^Уе — V. +__2_ дь^Л - (!)

г дг дг г2 де2 дг2 г2 г2 де

(д , , 1 д . , д . , 2

-Р — (Тет) + "ТГ" (Тее) + ТГ" (Тех) + "Тте

дг г дг дг г

Р I дЬ Ут дг г де

др

+ х дг ) дг + ^ \ г дг \ дг

1 д ( дь

1 д2ьг д 2ьх

+ -7Г + Х

г2 де2 дг2

1 д 1 д д

-Г - Р "7Т (гТхт) + "7Т (Тхе) + 7Г (Тхх Л ; г дг г де дг

д 1 д д

ТТ (гЬт) + -ТГ- (гЬе) + ТГ (гЬх) = 0, дг г де дг

где ьт, ье, — компоненты скоростии, р — давление, Ь — время, т^- — подсеточные

напряжения Рейнольдса, г, ] соответствуют г, е, г.

На рис. 1 отражена физическая область рассматриваемой задачи: К — радиус, Н —

высота, е = 2п — радиальный угол цилиндра. Лопасти пропеллера имеют форму прямо-

1 н2 — Н1 1 угольника, длина которого Ьрт = К2 — К1 = 4К, высота Нрх = —;-, ^ (Н2 + Н1)

2

вт а ' 2

(в

3 Н, время вращения за один период Тр, угол наклона лопасти а = агС^ ^ л ) , скорость

3

Л

пропеллера ьр. Лопасти пропеллера прикреплены к внутреннему цилиндру радиуса К1, а К2 — расстояние от центра внутри цилиндра до края лопасти, число лопастей К = 33. Используются следующие начальные (2) и граничные (3), (4) условия:

ьт (Ь = 0, г, е, г) = 0, ье (Ь = 0, г, е, г) = 0, (Ь = 0, г, е, г) = 0, ьр (Ь = 0, г, е, г) = 0.

Рис. 1. Физическая область: а — цилиндр; б — лопасть пропеллера

Уг

Уг

Уе Уг Уг

Уг

Уг

Ь, т = Е, е, г) = 0, у£ (Ь, т = Е, е,г) = 0, Уг (Ь, т = Е, е, г) = 0,

Ь, т = 0, е, г) = - (уг (Ь, Ат, е, г) + уг (Ь, Ат, е + п, г)), 2

Ь, т = 0, е, г) = у£ (Ь, т = Ат, е, г), Ь, т = 0, е, г) = Уг (Ь, т = Ат, е, г), Ь, т, е, г = 0) = 0, уг (Ь, т, е, г = Н) = 0, Ь, т, е, г = 0) = 0, Уг (Ь, т, е, г = Н) = 0, Ь, т, е, г = 0) = 0, Уг (Ь,т,е,г = Н) = 0.

(3)

При (Ех < т < Я2), (Н < г < Н2)

I — (ех < е < е2),

II — (е3 < е < е4),

III — (е5 < е < ее).

уг (Ь, т, е, г) = 0, у£ (Ь, т, е, г) = А, уг (Ь, т, е, г) = В,

ур (Ь, т, е, г) = у/(у22 (Ь, т, е, г) + у2 (Ь, т, е, г)), -

где А = 2,

(У.

в = -.

2

(4)

С целью устранения неопределенностей, появляющихся в основных уравнениях (1) при т = 0, произведено преобразование (рис. 2), позволяющее отобразить физическую область (т, е, г) на вычислительную (хх, х2, х3), имеющую форму куба [6], с использованием формул: т = ех1, е = Фх2, г = Нх3.

а

Рис. 2. Расположения лопастей физической области в вычислительной области в разные моменты времени вращения в сечении Хз = И\: а — Ь = ¿о; б — Ь = ¿о + ДЬ.

После проведения данного преобразования необходимо обезразмерить основные уравнения. Для этого выбираются характерные значения: тангенциальная скорость ^¿1; размеры дискретной области: ¿1 — длина области по оси ж1? ¿2 — длина области по оси ж2, ¿з — длина области по оси жз.

Используя указанные величины, введем следующие обозначения:

ж1

и1 =

ж1

¿1 ' и1

X о

т =у- , Р*

¿2 1

а (£1) = — ,

и2

_ £2 = «2

Р

'

X о

¿з'

«3 = 1 _ Йё =

V

'

(5)

где Йе — число Рейнольдса.

Подставляя безразмерные величины (5), получим:

дил _ дил ¿1_ (¿2^1) «3 д«1 2 + а ¿2^---+ а—«2^--1--7--77а--а ^¿1 («2)

дт дж1

= _а^Р + ¿* _!_

дж1 Йе

Ф

а

'дЖ2

д (дг71

¿з

+а2

д дж1

_!_ ¿1А

Н2 ¿3 дхз 1

дж1 \ дж1 д«1

+ а

Н джз 2 1 д

--Т11 + а

а

джз

¿1А

Ф дж

2

2 г 2 - 2 2 ¿1

_ а ¿1 «1 _ а

Ф2 ¿2 дж2 2 ¿1 д«2

Ф ¿2 дж2

д«1

дж2

+

+

( Т1 2) + ¿т^ (_Т1з) _ а ( Т22) '

Н джз

д«2 г _ д«2 ¿1_ д«2 (Ь2Ьl) из д«2

—--+ а ¿2«1---+ а—«2т,--1--7---

дт дж1 Ф дж2 ¿з

1 др 1

_а—---+ ¿2 —

Фдж2 2 Йе

д / ди

а

+

А ¿1

Н2 ¿3 джз

дж1

д«2 джз

дж

+ а

Н джз 1 ¿1 д

+ а ¿2^ («1 «2) д«2

Ф2 ¿2 дж2 V дж2

+

2г2 - , 2 2 ¿2 ди1

_а «2 + а — — -— 1 Ф ¿2 дж2

+

д ¿1 д ¿1 д

(_ Т21) + а--— ( Т22) + (_Т2з) + 2 а (_Т12)

дж

Ф дж2

Н джз

(6)

а

2

+

2

2

1

дия г _ дия дщ (¿2^1) и дщ

-дТ + и1 дХ1 + а ¥и2 дХ2 + НдХЗ

1 др

+Ь2— 2 Ие

д (ди,

а

дх1 V дх1

+ а2

1 д (ди

¥2 ¿2 дх2

дх2

+

Н дх3 1 д / ди

Н2 ¿3 дхз

дх3

+

+а2

д дх1

1

а

Т31

+ а —— ( Т3 2) + (-Т33) -

Ь дх2

Н дх3

Подставляя безразмерные величины (5) в уравнение неразрывности, получим

1 д (е- и,)+11 д

е2х1 ¿1дх

¥ еХ1 ¿2дх2

(и2) +

1 д Н ¿3дх3

(и3) = 0.

(8)

(9)

В (6)-(9) присутствуют подсеточные напряжения Рейнольдса, для моделирования которых используется вязкостная модель

_ 0 о

Т3--— Ткк = — ,

(10)

где турбулентная вязкость представляется в виде

^ = (СзА)2 (2&зБ

2 ~а \ 1/2 .

¿3) .

Б „ 2

'¿3

ш, + №

дх3 дх

6Ц

1,

— символ Кронекера; Cs — коэффициент, для данной задачи принят

г = ^

равным — 0.10 [4]; А = (еХ1 Ах1 Ах2 Ах3)1/3 — характерная длина фильтра, которая имеет порядок шага сетки; Ах1 = 1/^1, Ах2 = 1/^2, Ах3 = 1/^3 — шаг по вычисляемой сетке в направлении соответствующей оси (х1, х2, х3), N1, N2, N3 — количество узлов.

Таким образом, задача (1) заменена уравнениями (6)-(9), т.е. решение будет осуществляться в декартовой системе координат, при этом начальные и граничные условия для данной постановки соответствуют условиям (2)-(4). Следует учесть также выполнение условия периодичности для скорости пропеллера:

ур (т + Ат к, х1, х2 + Ах2, х3) = Ур (т, х1, х2, х3) , 0 < к < Тр,

которое свидетельствует об изменении расположения лопастей пропеллера при отображении их на декартову плоскость. При этом Ат = 1/Тр.

Компоненты скорости и давления, определенные из уравнений (6)-(9), используются для анализа тензора турбулентных характеристик. Турбулентные члены (и,и7) можно рассчитать лишь приближенно — на основании каких-либо гипотез. Опираясь на исследования [7], приведем результаты вычислении турбулентных членов:

(иН)

1 т " = (РТ)га+ 1 ,т,д 61 ,

п+ 2 , т , ч

61

-А

Ах21

дх1 / п+1 ,т,ч

'т

дх1 / п+1 ,т,ч

+

д^1

дх,

п+ 2 ,т,ч ^

ь

3

3

1

*

2

<и1и2>

2/ п,т+ 2 ,я = (Рт)га,т+ 1 82 г

>. Дх2

-А—2

-А 3

2

0£2

дх2 / . 1

2 / п,т+2

дЦ1 дх2

+ , ди

1 \ дхг

<иги3>,

Аги3>п,т,д+ 2

(Рт)

п,т+ 2

1 Оз — 'п,т,д+ 3 г

п,т+ 2 ,

ди

дхз

3 / п,т,д+ 2

дК

дх3

+

3 / п,т,д+ 2

дЦ3

дхг / п,т,д+ 2

(11)

где г = 1, 2, 3; 8^ = < ' — символ Кронекера, п — номер ячейки вдоль оси х1,

(О, г = ^

т — вдоль оси х2, д — вдоль оси х3; и1 = «1, и2 = «2, и3 = «3 — компоненты скорости, «г — пульсационные составляющие скорости, А = 0.02.

*

*

2. Численный метод

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

дробных шагов при использовании метода прогонки: и* - ип

1) —--= - (ипУ) и* - VДи*;

Дт

, Л Уи*

2) др = д7;;

ип+1 - и*

3) —ДГ~ =

На втором этапе по найденному промежуточному полю скорости определяется поле давления. Уравнение Пуассона для поля давления решается методом Фурье в сочетании с методом матричной прогонки, которая применяется для определения коэффициентов Фурье. На третьем этапе предполагается, что перенос осуществляется только за счет градиента давления [8-10].

3. Численные результаты

Численная модель позволяет описывать турбулентное движение нестационарных потоков в цилиндрической области. Расчеты были проведены для различных скоростей вращения лопастей внутри цилиндра. Для расчетов использовались сетки размером 20 х 20 х 20; 60 х 60 х 60. Шаг по времени, Дт, принимался равным 0.01; 0.001; 0.0001 соответственно. Указанными значениями Дт определяется характер движения жидкости: при Дт = 0.01 — ламинарное течение, Дт = 0.001 — переход от ламинарного течения к турбулентному, Дт = 0.0001 — развитое турбулентное течение.

На основе построенной модели определены следующие характеристики: турбулентная кинетическая энергия, корреляция, турбулентные масштабы, средние и пульсацион-ные характеристики течения для различных значений Рейнольдса, а также для разной численности лопастей пропеллера: Ке = 1000; 10 000, Дт = 0.01; 0.001, К = 3; 6; 12.

Установлено, что увеличение числа Яе приводит к появлению вихрей. Полученные результаты моделирования адекватно согласуются с экспериментальными данными, указанными в работах [1-5].

После получения расчетных данных для наглядности результатов произведено обратное преобразование:

1

еХх

1

1

-п

П = ф П

=

которое позволило отобразить на рис. 3 картину течения в виде изолиний скорости в сечениях £1 = 10°; = 200° при количестве лопастей К = 3 и угле наклона лопастей пропеллера а = 45°. Показан характер течения при А т = 0.001, т = 6000 Ат. Область белого цвета соответствует наибольшей скорости, черные области — наименьшей.

На рис. 4 и 5 представлено сравнение результатов расчета отношений радиальной пульсационной скорости п'г к скорости пропеллера (рис. 4) и кинетической энергии к квадрату скорости пропеллера (рис. 5) в зависимости от оси г, полученных в рамках данной задачи с экспериментом [4, 5]. Все расчеты проводились для ур = 0.09 при А т = 0.0001, на сечении, где значение г = 0.16, а угол £ = 150° в физической области. Поскольку расположение пропеллера, определенное для поставленной выше задачи (1), (3), не соответствует эксперименту, опубликованному в [4, 5], было произведено сопоставление расчетных данных перемещением (поворотом) физической области до совпадения с экспериментальной областью. Таким образом, осуществлено объективное сравнение эксперимента с полученными расчетными данными. Из рисунков видно, что отклонение расчетных данных от эксперимента небольшое, что свидетельствует об адекватности предложенного алгоритма.

пг =

Рис. 3. Изолинии в трехмерной области, при количестве лопастей К = 3 и угле наклона лопастей пропеллера а = 45°

Рис. 4. Отношение радиальной пульсацион-ной скорости к скорости пропеллера, ь!г. /гр, в зависимости от оси г: + — эксперимент [4, 5]; сплошная линия — результаты, полученные в настоящей работе

Рис. 5. Отношение кинетической энергии к квадрату скорости пропеллера: Е& — турбулентная кинетическая энергия; + — экспериментальные данные, штриховая линия — расчет ЬЕЯ(Я) [4]; сплошная линия — результаты, полученные в настоящей работе

На рис. 6 показаны рассчитанные по начальному полю коэффициенты корреляции скорости, КХ1Х1, КХ2 Х1, Яхз Х1 (в вычислительной области) (11), в различных точках в зависимости от расстояния между ними, п = 1, 2, N/2, N = N1 = N2 = N3, (при п > N/2). Характер изменения корреляций, изображенных на рис. 6, соответствует изменению коэффициентов корреляции, расчитанных в работе [7]:

К = 1

гХ1 2^ХЛ — п) ^3 (ц*.

'N-4 N N

XXX М ,ш,д

п=1 т=1 д=1

+

N N N

+ X/ X/ X/ (и^)п-7],т,д

п=1+п т=1 Ч=1

N N-n N

К1 х2

2^-1 (N2 — г) N3 (и2)

Ш Ш Ел ^

п,ш+п,д

+

п=1 т=1 д=1

МММ

+ X/ X [(и^)п,ш,д (и)п,т-7],д

п=1 т=1+п Ч=1

(12)

n n n-n

Кг Х3

2^-1 Nхх2 N3 — п) (и2)

Ш Ш Ш ^ п,т,д+ц

+

п=1 т=1 д=1

МММ

+ X/ X [(и^)п,ш,д (и^)п,т,д-7]

п=1 т=1 д=1+п

1

1

0.5 0.4 0.3 0.2 0.1 0.0 -0.1

Рис. 6. Коэффициенты корреляции: • — RX1Х1; ▲ — RX2 X1; ■ — Rx3 X1

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

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

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

Список литературы

[1] Derksen J.J., Van den Akker H.E.A. Large eddy simulations on the flow driven by a Rushton turbine // AIChE J. 1999. Vol. 45, N 2. P. 209-221.

[2] Yianneskis M., Popiolek Z., Whitelaw J.H. An experimental study of the steady and unsteady flow characteristics of stirred reactors // J. of Fluid Mechanics. 1987. N 175. P. 537-555.

[3] Piomelli U., Moin P., Ferziger J.H. Model consistency in large eddy simulation of turbulent channel flows // Physics of Fluids. 1988. N 31. P. 1884-1891.

[4] Hartmann H., Derksen J.J., Van den Akker H.E.A. Macro-instability uncovered in a Rushton turbine stirred tank // AIChE J. 2004. Vol. 50, N 10. P. 2383-2393.

[5] Verzioooy R., Iaccarino G., Fatica M., Orlandiz P. Flow in an impeller stirred tank using an immersed boundary method // AIChE J. 2004. Vol. 50, N 6. P. 1109-1118.

я я

_I___I___I_I_I_I_I_

0.0 0.1 0.2 0.3 0.4

[6] БЕЛОЦЕРКОВСКИй О.М. Прямое численное моделирование "переходных" течений газа и задач турбулентности // Механика турбулентных потоков. М.: Наука, 1980. C. 70-109.

[7] Иевлев В.М. Численное моделирование турбулентных течений. М.: Наука, 1990. 216 с.

[8] Флетчер К. Вычислительные методы в динамике жидкости. М.: Мир, 1991. 552 с.

[9] Яненко Н.Н. Метод дробных шагов решения многомерных задач математической физики. Новосибирск: Наука, СО, 1967. 196 с.

[10] Пейре Р., Тейлор Т.Д. Вычислительные методы в задачах механики жидкости. Л.: Гидрометеоиздат, 1986. 352 с.

Поступила в редакцию 26 февраля 2008 г., в переработанном виде —19 июня 2008 г.