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

УДК 621.928.37

МАТВИЕНКО ОЛЕГ ВИКТОРОВИЧ, докт. физ.-мат. наук, профессор, [email protected]

БАЗУЕВ ВИКТОР ПАВЛОВИЧ, ст. научный сотрудник, ЮЖАНОВА НАТАЛИЯ КОНСТАНТИНОВНА, аспирант, [email protected]

Томский государственный архитектурно-строительный университет, 634003, г. Томск, пл. Соляная, 2

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

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

Ключевые слова: битумные вяжущие; дорожное строительство; гидродинамика; закрученные потоки; турбулентность; реламинаризация.

OLEG V. MATVIENKO, Dr. Tech. Sc., Prof., [email protected]

VICTOR P. BAZUEV, senior researcher, [email protected] NATALIA K. YUZHANOVA, P.G., [email protected]

Tomsk State University of Architecture and Building, 2 Solyanaya sq., Tomsk, 634003, Russia

NUMERICAL INVESTIGATION OF TRANSITION TO THE TURBULENT REGIME OF INTERNAL SWIRL FLOWS OF BITUMINOUS BINDERS

The effect of swirl on the processes of turbulization and relaminarization of a bitumen binders stream in technological devices for manufacture of new materials for road construction has been investigated. It was determined that the occurrence of a reverse-flow zone prevents turbulent fluctuations around the axis due to a conservative nature of the effects of centrifugal force on the parameters of turbulence in this region of flow. Thereof, the laminar flow regime is maintained in the flow core.

Key words: bituminous binders; road construction; fluid dynamic theory; swirl flows; turbulization; relaminarization.

© О.В. Матвиенко, В.П. Базуев, Н.К. Южанова, 2013

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

В работах [1, 2] авторами исследован процесс приготовления нового вяжущего с применением принципа кавитационно-смесительного диспергирования, но в некоторых случаях необходимо применять другие технические решения, а именно переход к турбулентному режиму течения в закрученных потоках битумов или совместное смешение с применением кавитационно-смесительного диспергирования для получения качественных новых вяжущих. Это особо важно при приготовлении полимерных битумных вяжущих, т. к. полимеры обладают сложными структурными связями и технологически сложно получить качественное полимерно-битумное вяжущее.

Математическая модель

Одним из подходов исследования перехода к турбулентности является анализ распределения в потоке кинетической энергии турбулентности и ее зависимости от возмущений, возникающих в потоке. В рамках этого подхода для описания поля течения в цилиндрическом канале используются двумерные осе-симметричные уравнения Рейнольдса, записанные относительно осредненных по времени составляющих скорости: и , V, ^ и давления р [3, 4]:

дри + 1 дршг дх г дх

■ = 0.

(1)

дри2 1 дришг

дх

дг

др дх

д_

дх

€

( 2 ди 2 ( ди

+1—

г дг

дриш + 1 дрш2 г дх г дг

1 д

+--

г дг

( ди

* [тг

дх 3[ дх

дш дх

1 дшг г дг

(2)

др дг

д_

дх

дш дх

ди дг

€

„дш 2 (ди 1 дшг

2---^ +--

дг 3 [ дх г дг

дрим> дх

1 дршлг г дх

д дх

дw

тх

1. д

г2 дг

_ 2

г 3 Та ^

дг [ г

(3)

(4)

В отсутствии турбулентных пульсаций осредненные по Рейнольдсу параметры совпадают с мгновенными, и уравнения Рейнольдса переходят в уравнения Навье - Стокса для описания ламинарного течения.

Исследование характеристик турбулентности осуществлялось с использованием двупараметрической к-г модели, адаптированной Джонсом и Лаундером для расчета течений с низкими числами Рейнольдса [5-7]:

дрик 1 фукг

дх

дх

д ц е/ дк' 1 д +

дх _Ск дх г дг _Ск дг _

дрие 1 друвг д дх

дх

дх

ц / де 1 д де"

_Се дх_ г дг _ Се дг

■Ок-рг-В, (5)

+ (С1 - С2рг)~ + Е . (6)

В формулах (5) и (6) Ок - диссипативная функция, которая рассчитывается как

Ок = ц I2

В = 2ц

V ^ +(ди + ду г) \дг дх

(д4кI2 (д4к

дw дх

+ 1 г

дw / г дг

Е = 2

(д 2и >

дх2

дх

ди ^

J \ 2

I г— дг I дг

дх

(д ^ >

(7)

(8)

дх2

1 (д_( _дv^ дг I дг

"(д2 ^ ^

дх2

1 ( д_( дг I дг

(9)

Значения констант выбираются в соответствии с рекомендациями [8]: С = 1,44 (1 + С4Ш е ) , С2 = 1,92 (1 - СзШ,), Сц= 0,09, Ск = 1, ^ = 1,3. Турбулентное и градиентное числа Ричардсона определяются следующим выражением:

д (

к2 w д(wг)

Я =■

g

Ок дг V г

(10)

е" г~ дг

Эффективная вязкость (це/) определяется как сумма молекулярной (ц) и турбулентной вязкости (ц). Турбулентная вязкость может быть рассчитана с использованием к-г модели турбулентности [9]:

Ц = Сц /црк 2е-1, (11)

где Сц = 0,09 - константа модели турбулентности,

/ц = ехр (-3,4/(1 + 0,02Яе )2) , /2 = 1 - 0,3ехр (-Яе^ ) - функции модели турбулентности, зависящие от турбулентного числа Рейнольдса = рк2 /це .

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

Граничные условия на входе определяются для всех переменных. Задается распределение скорости потока, а кинетическая энергия турбулентности на входе берется пропорциональной кинетической энергии осреднен-ного течения.

, 2г '

х = 0 : и = ит, V = 0 , м = и.О.

к = Ти (и2 + w_), е = 2£3/2/М).

\ 1п 1п у ' т VI/

Здесь ^ = 0,005 ; Ти - константы модели; ё - диаметр канала; О, п - параметры, характеризующие интенсивность закрутки потока на входе в канал.

При п = 1 реализуется закрутка по закону вращения твердого тела, и параметр О характеризует величину угловой скорости а = 2Оит/ё. При

п = 0 моделируется закрутка с постоянным углом ф = аг^ (3/2 О) .

Влияние закрутки на структуру течения и теплообмен в потоке удобно характеризовать интегральным параметром закрутки Хигира - Бэра [9], представляющим собой отношение осевой составляющей потока момента количества движения к произведению радиуса канала и осевой составляющей потока количества движения:

12 К ё%2 )

Ф. = I рuwг2ёг — I ри2гёг . (12)

0 / [ 2 0 )

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

ф.,^- (13)

п + 3

В настоящей работе будем использовать Ф. т для сравнения потоков,

закрученных с использованием различных законов закрутки.

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

т : ди 0 0 дм 0 дк 0 де 0

х = Ь: — = 0, ш = 0, — = 0, — = 0, — = 0.

дх дх дх дх

На оси канала задаются традиционные условия:

„ = 0: * = 0, ш = 0, м = 0, д± = 0, * = 0. дг дг дг

На стенках канала выполняется условие прилипания. Для определения турбулентных характеристик предполагается локальное равновесие в пристеночной области:

г = ё/ 2 : и = 0, ш = 0, м = 0, к = 0, е = 0.

Численная реализация

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

При аппроксимации конвективных членов применялись разности против потока с использованием схемы [10]. При моделировании диффузионных членов использовалась центрально-разностная аппроксимация.

Система конечно-разностных уравнений является нелинейной, и для ее решения применялся итерационный метод. На каждой итерации использовалась продольно-поперечная прогонка. Давление рассчитывалось с помощью итерационных процедур 8ТМРЬБС [11].

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

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

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

Г ( 2"

и (х = 0, г ) = 2ип

1 -'2г

(14)

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

Возмущения при Яе < Яе* могут не только вноситься в канал из области х < 0, но и при определенных условиях втекания в канал генерироваться на входе. Рассмотрим режим гидродинамически развивающегося течения, когда на входе в канал задан плоский профиль осевой скорости со скачком на стенке:

и (х = 0, г < й/2) = ип, и (х = 0, г = й/2) = 0 .

Такой режим течения возможен при напорном истечении из емкости большого диаметра В в трубу диаметром й << В . Наличие острой кромки в точке х = 0, г = й/2 приводит к разрыву скорости и возникновению возмущений. Последние, возникнув в окрестности острой кромки, сносятся

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

0,10 0,08 0,06 0,04 0,02

Осевая координата х, м

Рис. 1. Изолинии турбулентной кинетической энергии. иП = 0,1 м/с, Ти = 0,5 , Ке = 1360

При иП > и, =ц0/pd ■ Ке, турбулентность в потоке обусловливается уже не только первоначальной турбулизацией потока или возмущениями на входе, но и передачей энергии от осредненного течения турбулентному через напряжения Рейнольдса. В результате процесс реламинаризации отсутствует и изолинии турбулентной кинетической энергии перестают быть замкнутыми. Возникает так называемый режим пристенной турбулентности (рис. 2).

Осеваякоордината х, м

Рис. 2. Изолинии турбулентной кинетической энергии. и1п = 0,2 м/с, Ти = 0,5 , Ке = 2720

С дальнейшим увеличением иП область, занимаемая турбулентным течением, становится шире и распространяется в ядро потока, постепенно занимая всю область течения. Наступает режим развитой турбулентности.

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

ющая ко входу в канал. Ее наличие объясняется подачей в канал потока, возмущенного в области х < 0. Интенсивность этого возмущения можно характеризовать параметром к т = Ти • и2 . Расчеты показывают, что с увеличением кт область, занимаемая возмущенным течением, продвигается вниз по потоку. Однако область первоначально турбулизированного течения является ограниченной некоторым предельным значением х1, за которым происходит реламинаризация течения.

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

Влияние закрутки на турбулизацию течения

Закрутка потока приводит к появлению в потоке тангенциальной составляющей скорости м и формированию центробежной силы ** = рм2 /г .

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

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

Влияние закрутки на распределение кинетической энергии турбулентности иллюстрирует рис. 4.

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

й Рч

.08

б

,06

,04

,02

0,2

Осевая координата х, м

Осевая координата х, 1м

Осевая координата х, м

Рис. 3. Линии тока ип = 0,1 м/с, Ке = 1360 , Ти = 0,03, Ф. = 8,5 : а - п = 0, б - п = 0,5, в -п = 1, г- п = 2

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

0,10 -0,08 -0,06 -0,04 -0,02 -

0,4 0,6

Осевая координата х, м

й

н 0,06 й

В 0,040

I 0,02

О О

к

0,4 0,6 Осевая координата х, м

Осевая координата х, м

Р ^ о и

3 0,10 -

с 0,08 -

а

нат 0,06 -

дин 0,04 -

рд 0,02 -

о

0,001 —0,004

0,4 0,6

Осевая координата х, м

0,4 0,6

Осевая координата х, м

~ 0,08 £

Й 0,06

0,04 0,02

0,4 0,6

Осевая координата х, м

Рис. 4. Изолинии турбулентной кинетической энергии ип = 0,1 м/с; Ке = 1360; Ти = 0,03; а - п = 0 ; Ф, = 1,25 ; б - п = 0, Ф, = 8,5; в - п = 1, Ф, = 1,25 ; г -п = 1, Ф, = 8,5 ; д - п = 2 , Ф, = 1,25 ; е - п = 2 , Ф, = 8,5

а

0,2

0,8

0,10

б

0,08

0,2

0,8

в

г

0,2

0,8

д

0,2

0,8

е

0,8

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

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

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

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

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

Около стенки турбулентное число Ричардсона принимает отрицательные значения, что свидетельствует о турбулизации течения. По мере удаления от входа наблюдается распространение зоны активного воздействия поля центробежной силы из пристеночной области в ядро потока. С увеличением интенсивности закрутки происходит рост как мощности центробежной силы, так и турбулентности. При слабой и средней закрутке Ф, < 1 прирост мощности турбулентных возмущений превышает увеличение мощности центробежной силы, что приводит к увеличению зоны отрицательных значений турбулентного числа Ричардсона и увеличению с закруткой максимального абсолютного значения К. Однако с появлением рециркуляции мощность центробежной

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

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

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

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

Библиографический список

1. Базуев, В.П. Моделирование процесса модифицирования битума в кавитационно-смесительном диспергаторе / В.П. Базуев, О.В. Матвиенко, В. Л. Вороненко // Вестник ТГАСУ. - 2010. - № 4. - С. 121-128.

2. Матвиенко, О.В. Исследование динамики пузырька в закрученном потоке нелинейно-вязкой жидкости / О.В. Матвиенко, М.В. Агафонцева, В.П. Базуев // Вестник ТГАСУ. - 2012. -№ 4. - С. 144-156.

3. Матвиенко, О.В. Численное исследование сепарационных характеристик гидроциклона при различных режимах загрузки твердой фазы / О.В. Матвиенко, И.Г. Дик // Теоретические основы химической технологии. - 2006. - Т. 40. - № 2. - С. 219-224.

4. Дик, И.Г. Некоторые закономерности теплообмена внутренних закрученных потоков / И.Г. Дик, О.В. Матвиенко // Известия СО АН СССР. Сер. техн. наук. - 1989. - Вып. 3. -С. 40-43.

5. Jones, W.P. The calculation of low Reynolds number phenomena with a two-equation model of turbulence / W.P. Jones, B.E. Launder // Int. J. of Heat Mass Transfer. - 1973. - No 16. -P. 1119-1130.

6. Матвиенко, О.В. Анализ моделей турбулентности и исследование структуры течения в гидроциклоне / О.В. Матвиенко // Инженерно-физический журнал. - 2004. - Т. 77. -№ 2. - С. 58-64.

7. Wilcox, D.C. Formulation of the k -ю Turbulence Model Revisited // AIAA Journal. -2008. - No 46 (11). - P. 2823-2838.

8. Piquet, J. Turbulent Flows: Models and Physics / J. Piquet. - Berlin : Springer, 1999.

9. Гупта, А. Закрученные потоки / А. Гупта, Д. Лилли, Н. Сайред. - М. : Мир, 1987.

10. Ferziger, J.H. Computational Methods for Fluid Dynamics / J.H. Ferziger, M. Peric. - Berlin : Springer, 1996.

11. Versteeg, H.K. An Introduction to Computational Fluid Dynamics / H.K. Versteeg, W. Mala-lasekera. Longman. - London, 1998.

REFERENCES

1. Bazuev, V.P., Matvienko, O.V., Voronenko, V.L. Modelirovanie processa modificirovanija bituma v kavitacionno-smesitel'nom dispergatore [Modeling of the process of modifying the bi-

tumen in cavitation-mixing disperser] // Vestnik of Tomsk State University of Architecture and Building. - 2010. - No. 4. - P. 121-128. (rus)

2. Matvienko, O.V., Agafontseva, M.V., Bazuev, V.P. Issledovanie dinamiki puzyr'ka v zakru-chennom potoke nelinejno-vjazkoj zhidkosti [Investigation of dynamics of bubble in the swirling flow of nonlinear viscous fluid] // Vestnik of Tomsk State University of Architecture and Building. - 2012. - No.4. - P. 144-156. (rus)

3. Matvienko, O.V., Dik, I.G. Chislennoe issledovanie separacionnyh harakteristik gidrociklona pri razlichnyh rezhimah zagruzki tverdoj fazy [Numerical study of the hydrocyclone separation characteristics at different modes of solids charge] // Teoreticheskie osnovy himicheskoj tehnologii [Theoretical Bases of Chemical Engineering]. - 2006. - V. 40. - No. 2. - P. 219-224. (rus)

4. Dik, I.G., Matvienko, O.V. Nekotorye zakonomernosti teploobmena vnutrennih zakruchennyh potokov [Some regularities of heat transfer of internal swirling flow] // Izvestija SO AN SSSR [Bulletin ofAcademy of Sciences of the USSR]. Ser. tehn. nauk. - 1989. - No. 3. - P. 40-43. (rus)

5. Jones, W.P., Launder, B.E. The calculation of low Reynolds number phenomena with a two-equation model of turbulence // Int. J. of Heat Mass Transfer. - 1973. - No. 16. - P. 1119-1130.

6. Matvienko, O.V. Analiz modelej turbulentnosti i issledovanie struktury techenija v gidroci-klone [Analysis of turbulence models and research of structure of flow in a hydrocyclone] // Inzhenerno-fizicheskij zhurnal [Journal of Engineering Physics]. - 2004. - V. 77. - No. 2. -P. 58-64. (rus)

7. Wilcox, D.C. Formulation of the k -m Turbulence Model Revisited // AIAA Journal. -2008. - No. 46 (11). - P. 2823-2838.

8. Piquet, J. Turbulent Flows: Models and Physics. - Berlin, Springer, 1999.

9. Gupta, A., Lilli, D., Sajred, N. Zakruchennye potoki [Swirling flows]. - Moscow, Mir, 1987. (rus)

10. Ferziger, J.H., Peric, M. Computational Methods for Fluid Dynamics. - Berlin, Springer, 1996.

11. Versteeg, H.K., Malalasekera, W. An Introduction to Computational Fluid Dynamics. Longman. - London, 1998.