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

2016 Математика и механика № 6(44)

УДК 532.5

Б01 10.17223/19988621/44/8

В.В. Чуруксаева, А.В. Старченко

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

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

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

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

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

1 Работа выполнена при поддержке РФФИ и Администрации Томской области в рамках научного проекта № 16-41-700178 р_а.

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

Рассматривается двухфазное изотермическое движение смеси «вода - легкие частицы» в открытом канале или русле реки. Межфазовый обмен массой и теплом не учитывается в силу близких значений температуры воды и окружающей среды и их незначительных изменений за период моделирования. Плотность льда р0 = 910 кг/м3 меньше плотности воды р0 = 1000 кг/м3, и потому считается, что ледяные частицы плотно расположены в приповерхностном слое воды и их концентрация остается постоянной на входе в канал (или рассматриваемый участок реки).

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

Рис. 1. Физическая постановка задачи Fig. 1. Physical statement of the problem

Для математического описания данного процесса будем использовать уравнения механики взаимодействующих взаимопроникающих континуумов [1]. Запишем уравнения, описывающие движение жидкой фазы:

Ф i + ФЛ dt

dx,

= 0;

(1)

dPiw ij Ф i wikwi,

- +-— = -a,

dr+pi gj kj +Ttkj)]+Su;

dt

dx.

dx.

dx.

] = 1,2,3. (2)

По повторяющемуся индексу (к=1, 2, 3) проводится суммирование. Здесь индекс «/» относится к фазе воды; х1, х2, х3 - координаты декартовой системы;

g = ( 0,0, —g) - вектор ускорения свободного падения; wi =(Ч

И' "72' i3

) - вектор

скорости движения воды; / - время; - источниковый член, описывающий силы взаимодействия фаз и влияние силы Кориолиса.

Z

X

р, = р°а,, где р° - истинная плотность воды (принимается равной

1000 кг/м3), а, - объемная доля воды 0 <а, < 1.

тк = ЦI

дХ:

V ]

дхк

ний несущей среде,

(

г _ г

т1к: = цI

^к + .5Р

Vдx]

дхк

div Р, - компоненты тензора вязких напряже-

- з р{д к]ке; - компоненты тензора турбулентных напря-

жений,

- молекулярная вязкость воды; Ц - турбулентная вязкость воды; 8к]- -символ Кронекера; ке - кинетическая энергия турбулентности.

Уравнения для дисперсной фазы (континуума ледяных частиц) будут иметь следующий вид:

дРг , дРгРгк

дг дхк

■ = 0;

(3)

дРгРг]_ + дРгРгкРг] =-а др

дг

дх.

дхр + [аг ( + Тк )] + % ;

дх.

(] = 1, 2, 3).

(4)

По повторяющемуся индексу проводится суммирование (к=1, 2, 3). Здесь ин-

декс «г» относится к дисперсной фазе ледяных частиц; т

т*

г] гк]

- компоненты тен-

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

(

тгк] = Ц

V дх] дхк ,

2

- 3 ] div Рг . Рг = Р.Ч

где р° - истинная плотность льда ( р0 = 910 кг/м3), аг - объемная доля фазы ледяных частиц 0 <аг < 1, аг + а, = 1.

Рассмотрим источниковые члены, входящие в уравнения движения фаз. Источник £г- = + Ёум + Ёс в уравнении движения дисперсной фазы

представляет собой сумму следующих сил [1]: - сила Архимеда

РА = Р° аг I -ё I = РА1 +РА2 ; РА1 = Р аг ; РА2 = агё ,

щ

щ

вг

др дг

- + Р

дР дх.

I .

- сила вязкого трения, действующая на частицы со стороны несущей среды [2],

3 0 -1.65

Р = 3 ^ща^ - р | ( - р), а, > 0.8,

4 dt.fi

р» =

150

а2 (1 -а, +175 а р0 -

Ш )2

( - $) а, <08

где св = тах

24

Яе Р

(1 + 0.15Л°'687 ),0.44

- безразмерный коэффициент сопротив-

. . ^ \м>, - #,|^р0а, ления; а. - эффективный диаметр ледяных частиц, Яе = —-0— ;

- сила присоединенной массы, возникающая благодаря инерции в несущей среде,

0

РУМ = СУЫа грг

СуМ и 0.5 для сферических частиц;

0 |

вг вг

0 >

Ра - Р = -ЯР? (X, + (£0 -1 (1 -а,)

1Р° V

- сила Кориолиса.

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

(Н + Ч - хз )

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

Пусть а,, а.. - массовые доли плавающих ледяных частиц и воды соответственно. Тогда 0 <а, < 1, 0 <аг < 1, а, + аг- = 1;

Н+гь

| агаг = а1Д. = Н";

Н-Н.+1ь

Н+1ь Н+1ь

Н' = | агсЪ = а¡Н = | (1 -а.. )Ог = Н - Н".

Таким образом, проинтегрировав соотношение а, + а.. = 1 по глубине потока

1 ч+Н

Н, получим а¡Н + аД. = Н , где а, = — | а,аХ3. Здесь Н. - эффективная глубина

ч

слоя ледяных частиц, Н, - глубина слоя несущей фазы под слоем льда.

Предполагая, что распределение динамических параметров двухфазного течения по глубине близко к однородному, будем интегрировать записанные уравне-

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

дк" дк'% дк'%

+

дг дх

ду

= 0;

(5)

дк'% дк" й^ дк'%Л - +-- + -

= -^0 ^

Го д -1

vРz

дг дх

(1 -а,)

ду

ду

дУ дх

д(ъ + к )

дх

2 д

3 дх

2 К +^г )к "

дх дйг д% дх ду

диг'

дх _ +

1Чр° , „(ди, _ди1 _ ди1 ,

+ Р0 к (и, - й ) + (С^ + 1)-0 к ^ + и, + V, -¡У |-

р

Рг . „ (диг _ диг _ диг _

-Сум^к I ^ + иг + Vг | + /к " %г ,

дг

дх

ду

(6)

дк'% дк"и%, дк'%2

Р? Як"

( ~0 Л -1

чР,

дг

(1 -а,)

ду

ду

2 (г + ^гг )к"

ду ]

2д

3 ду д%,

дх

д(гъ + к) +

ду дх

<*+*' 'к' II + £ |]

V г к +2.1+к к

дх ду |

Р к" _

р° к" ( д% _ д% _ и%,

+ (еум +1)^4-I — + и,— + % —-1-с

дг

--Г V, -

дх ду

ум 0 Рг

дг

дх

ду

(7)

р, , „( д% _ д% _ д% , _ ' к I —^ + и—г- + %—г- I к - /к иг

и несущей фазы - воды

дк' дк'и, дк'%, — +-+-= 0;

дг дх ду

(8)

= - Як'

д +—

ду

а, +

0

^ -1

V Р/

дк'и1 дк'и, дк'и%1 дг дх

(1 -а,)

ду

д(ъ + к) +

дх дх

2 (V, +v,í )к'

дх

Р,

к ' (^ + ^ 1 + к к

, _ ч,, (ди, д%, ) 2 д

(V, + v,t )к' I—+ — |---

V ду дх 3 дх _ V дх ду I Р , „!- - \ , „( ди, _ ди, _ Щ }

+ к" (иг - щ ) - Сумк" + и + % -у | +

.„, диг _ диг _ диг . ,_ ,_ ._

+Смк !"дГ+и "дх+% "ду 1+'/к % -С/ ;

дх

дк'у дИугиг дН'у'2

дг

дх

ду

= - Я*

а■ +

(р? Р?

-1

(1 -а г)

д(ъ + к ) +

ду дх

, _ диг

^ дх ду )_

ду

2 (VI + Уй )А'дУ

дУ _

2 _д_

з ду

v¡h,(дUL + ^ 1 + кк

\ дх ду )

(10)

+к (V - V)- оумЬ"\^г + и■ + V | +

Рг

дг

дх

ду

, „(ду _ ду; )

+симк ^+и; +"дУ7^ - /И'и1 - С/1^ IV'

Здесь А - глубина всего потока, А"(г, х, у) - глубина слоя дисперсной фазы, (г, х, у), у; (г, х, у) - осредненные по глубине значения компонент вектора горизонтальной скорости м>; =(и;, ); 2ъ (х, у) - рельеф дна; р0, р? - истинные плотности льда и воды соответственно, я = 9.81м / с2 - ускорение свободного падения; ki - осредненная по глубине кинетическая энергия в слое дисперсной фазы; v;, уй - молекулярная и турбулентная вязкость дисперсной фазы; / - параметр

Кориолиса; С/ =

ЯП

И

0.333

- трение жидкости о дно канала (реки), п > 0 - коэффици-

ент Маннинга.

Для расчета турбулентных характеристик двухфазного течения используется высокорейнольдсовая к -е -модель турбулентности для осредненных уравнений [3], с модификацией Поурахмади и Хумфри [2] для учета влияния дисперсных частиц.

Начальные и граничные условия

В начальный момент времени г = 0 используются следующие значения параметров течения:

к" = к;се - известная величина;

о> ;

к = к - кс

иг~и1 о> о-

На входе в расчетную область значения параметров жидкой фазы и фазы частиц считаются известными, на выходной границе используется равенство нулю производных по внешней нормали к границе. Кроме того, на границах потока с берегом рассматривается трение как для жидкости, так и для частиц. На боковых стенках русла применяются условия непротекания и прилипания для компонент скорости. В случае, когда вблизи боковых стенок в несущей фазе справедливо соотношение vlt » vl, трение и турбулентные характеристики в несущей фазе в

пристеночной области определяются с помощью метода пристеночных функций

Лаундера - Сполдинга [5]. Трение фазы частиц о дно реки в прибрежной зоне и на

-1

участках отмелей выражается зависимостью с у , где с у = 0.0025 | [6].

Численный метод решения уравнений модели

Уравнения модели дискретизируются с использованием метода конечного объема [7] на разнесенных сетках (рис. 2), т. е. конечные объемы для компонент скорости сдвигаются на полшага сетки относительно центра Р конечного объема, используемого при аппроксимации уравнений для скалярных величин

И', И', к, е .

Рис. 2. Сеточный шаблон разностной схемы. Большими буквами отмечено положение центров конечных объемов, малыми - середин их граней [7]

Fig. 2. Mesh pattern of the difference scheme. Uppercase and lowercase letters indicate the centers of finite volumes and midpoints of their edges, respectively [7]

Конвективные слагаемые уравнений аппроксимируются с применением монотонных разностных схем высокого порядка (MLU [8] или MUSCL [9]). Диффузионные слагаемые представляются на основе центрально-разностной схемы второго порядка, при этом для снижения существенного ограничения на шаг интегрирования по времени для уравнений движения при аппроксимации членов, отвечающих за динамическое взаимодействие фаз (сила трения), применялась неявная аппроксимация и специальная процедура решения сеточных уравнений, позволяющая использовать построенную разностную схему для областей двухфазного течения, где отсутствует дисперсная фаза (к" = 0).

Кратко опишем предлагаемый подход.

Рассмотрим дискретный аналог уравнений (6) и (9) для внутреннего узла e расчетной сетки, представленной на рис. 2:

=ф0 +Рк;0 ( -Sie); (11)

Y 0 +4 PkJ 0 (ie - Se ). (12)

he - he

At

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

В правых частях сеточных уравнений (11), (12) отдельно выделены слагаемые

о

ph0°(íj0 -íle) и —0ph"0 (ít¡e -íi0), описывающие динамическое взаимодействие Р.-

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

0.5 Ах Ау

вие Куранта т <-—-—-¡=-) ограничение на шаг по

max|í0| Ау + max|vn| Ах + уgh (Ах + Ау)

времени при проведении расчетов течений с частицами, обладающими высокой динамической инерционностью. Кроме того, в узлах сетки, где h" = 0, система уравнений (11)-(12) имеет особенность (уравнение (12) обращается в тождество). В связи с этим перепишем систему (11), (12) в следующем виде:

(h'e 0 + Аф/г0' 0 )ще -А/Щ %е = А/Ф0 +h'e ; (13)

А „0 Л „0

1 + £0 Аф

. Pi

íie -Р0Афи^ = А/Т0/К° + Í0. (14)

Pi

В (14) в ТИ'0 при И'0 <е, где е>0 - малая положительная величина, Т? принимается равной нулю. Определитель системы отличен от нуля и система решается по правилу Крамера. Аналогичным способом представляется система для компонент скоростей фаз у1п, У;П .

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

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

Результаты численных расчетов

Течение в и-образном канале

Экспериментальная установка представляет собой И-образный канал с дном в виде лотка со скошенными стенками (рис. 3). Экспериментальные исследования течения в данной установке представлены в работах [12,13]. В данном случае рассматривается эксперимент, описанный в [12]. Глубина потока на входе в канал

равнялась 0.45 м, расход воды на входе в канал - 0.16 м3/с . Так как данные о шероховатости стенок не представлены, коэффициент Маннинга выбирался равным 0.01413, что соответствует слабошероховатым стенкам (из бетона или пластика). Данные параметры соответствуют турбулентному течению с числом Рейнольдса, вычисленным по глубине потока Яей = 77170, Ег = 0,081. Для имитации ледяных частиц использовались сферические бусины из полипропилена диаметром ^,=0.005 м. Плотность материала частиц р0 = 900 кг/м3 близка к плотности льда. При проведении физического моделирования движения плавучих частиц они сбрасывались в стационарный поток с одинаковым расходом в конце прямого предвключенного участка (отмечен горизонтальной стрелкой на рисунке 3а). Измерения толщины слоя частиц и глубины потока в эксперименте проводились с точностью до 1 мм.

a b

Рис. 3. Канал с разворотом течения (размеры указаны в метрах): (а) - вид сверху; (b) - поперечное сечение Fig. 3. U-shaped laboratory flume (dimensions are in meters): (a), plan view, and (b), cross section

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

В результате расчетов обнаружено, что течение в канале с плавным разворотом с легкими частицами имеет сходные характеристики с течением в реке во время ледохода [12].

URROZ AND ЕТТЕМА

Рис. 4. Профили фронта слоя частиц по мере движения по каналу

(a) - эксперимент [12], (b) - расчет Fig. 4. Front profiles of the layer of particles while moving in the flume: (a) experimental [12] and (b) calculated data

Рис. 5. Мгновенные профили глубины h" (a) и модуля скорости дисперсной фазы (b) Fig. 5. Instantaneous profiles of the (a) depth h" and (b) velocity module of the dispersed phase

Течение в канале с резким поворотом

Рассматривается движение смеси «вода - ледяные частицы» в канале с резким поворотом. Геометрия канала изображена на рис. 6. Подробное исследование структуры однофазного течения в нем проводится в [3,10].

0.86 I "

ш in in

Резкое снижение уровня дна

¿2

4.43

Рис. 6. Лабораторная установка - канал с поворотом под углом 90° Fig. 6. Laboratory facility, a flume with a 90° bend

Входной участок канала имеет длину 5.555 м, ширину 0.86 м и ровное дно. Сразу перед поворотом уровень дна понижается на 0.013 м. Выходной участок имеет длину 4.43 м, ширину 0.72 м и ровное дно. На рисунке также отмечены области, где образуется возвратное течение, выявленные в процессе исследования однофазного течения в данном канале [3].

Рассмотрим течение воды с ледяными частицами со скоростью на входе в канал равной U0 = 0.2 м/с, начальной глубиной потока равной h = 0.175 м. Глубина слоя дисперсной фазы h" = 0.04 м. Параметры дисперсной фазы для базового и методических расчетов приведены в табл. 1.

Таблица 1

Параметры расчетов: 1 - базовый расчет, 2, 3 - методические расчеты

№ расчета Диаметр частиц, d , м Коэффициент формы, f Вязкость в слое частиц, v i, м2/с

1 0.1 0.166 0.01

2 0.01 0.166 0.01

3 0.1 1 0.01

Оценка влияния размера частиц

На рис. 7 представлено сравнение модуля скорости несущей фазы (а), модуля скорости дисперсной фазы (Ь), глубины слоя дисперсной фазы (с) для расчетных случаев (1) и (2) (см. табл. 1).

3

4

5

3

4

5

3

4

5

Рис. 7. Поля скорости несущей фазы (a) скорости дисперсной фазы (b) и глубины слоя

дисперсной фазы (с) для базового варианта (1) и расчета с мелкими частицами (2)

Fig. 7. Fields of the (a) carrier phase velocity, (b) dispersed phase velocity, and (c) depth of dispersed phase layer: (1), basic variant, and (2), computation with small particles

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

1. Скорости движения мелких и более крупных частиц существенно не различаются в области поворота течения. За поворотом у правой по течению стенки наблюдается зона рециркуляционного течения, более интенсивного для частиц с диаметром d = 0.1 м.

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

2. Более заметное влияние на скорости движения фаз и толщину слоя частиц оказывает небольшое изменение рельефа дна. Из рисунков видно, что изолинии модуля скорости жидкости и фазы мелких частиц резко изменяются вблизи снижения рельефа дна, т.е. малоинерционные частицы так же, как и жидкость, быст-

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

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

Оценка влияния коэффициента формы частиц

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

Рис. 8. Поля скорости несущей фазы (a) скорости дисперсной фазы (b) и свободной поверхности (с) для базового варианта (1) и расчета со сферическими частицами (3) Fig. 8. Velocity fields of the (a) carrier phase, (b) dispersed phase, and (c) free surface: (1), basic variant, and (3), computation with spherical particles

3

4

5

3

4

5

3

4

В области поворота течения наблюдается повышение уровня свободной поверхности в углу канала (более значительное для сферических частиц f = 1 ) и снижение за поворотом в области рециркуляционного течения. Такой характер распределения свободной поверхности в области поворота течения соответствует известным экспериментальным данным [14].

Выводы

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

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

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

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

ЛИТЕРАТУРА

1. Нигматулин Р.И. Динамика многофазных сред. Ч.1. Москва: Наука, 1987. 464 с.

2. Бубенчиков А.М., Старченко А.В. Численные модели динамики и горения аэродисперсных смесей в каналах. Томск: Изд-во Том. ун-та, 1998. 236 с.

3. Чуруксаева В.В., Старченко А.В. Математическая модель и численный метод для расчета турбулентного течения в русле реки // Вестн. Том. гос. ун-та. Математика и механика. 2015. № 6(38). С. 100-114. DOI: 10.17223/19988621/38/12.

4. Роди В. Модели турбулентности окружающей среды // Методы расчета турбулентных течений. М.: Мир, 1984. С. 276-278.

5. Launder B.E., Spalding D.B. The numerical computation of turbulent flows // Computer Methods in Applied Mechanics and Engineering. 1974. V. 2. No. 3. P. 269-289. DOI: 10.1016/0045-7825(74)90029-2.

6. Gidaspov D. Multiphase Flow and Fluidization: Continuum and Kinetic Theory Descriptions. Boston: Academic Press, 1994. 457 p.

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

8. van Leer B. Towards the Ultimate Conservative Difference Scheme, V. A Second Order Sequel to Godunov's Method // Journal of Computational Physics. 1979. No. 32. P. 101-136. DOI: 10.1016/0021-9991(79)90145-1.

9. Cada M., Torrilhon M. Compact third-order limiter functions for finite volume methods // J. Computational Physics. 2009. V. 228. P. 4118-4145. DOI: 10.1016/j.jcp.2009. 02.020.

10. Cea L., Puertas J., and Vazquez-Cendon M.E. Depth averaged modelling of turbulent shallow water flow with wet-dry fronts // Archives of computational methods in engineering. September 2007. V. 14. No. 3. P. 303-341. DOI: 10.1007/s11831-007-9009-3.

11. Hou J., Simons F., Mahgoub M., and Hinkelmann R. A robust well-balanced model on unstructured grids for shallow water flows with wetting and drying over complex topography // Computer methods in applied mechanics and engineering. 2013. No. 257. P. 126-149. DOI: 10.1016/j.cma.2013.01.015.

12. Urroz G.E., Ettema R. Bend ice jams: laboratory observations // Canadian Journal of Civil Engineering. 1992. V. 19. P. 855-864. DOI: 10.1139/l92-097.

13. Tsai W.F., Ettema R. Ice cover influence on transverse bed slopes in a curved alluvial channel // J. Hydraulic Research. 1994. V. 32. No. 4. P. 561-581. DOI: 10.1080/00221686. 1994.9728355.

14. Han S.S., Ramamurthy AS., and Biron P.M. Characteristics of Flow around Open Channel 90° Bends with Vanes // Journal of Irrigation and Drainage Engineering. 2011. V. 137. No. 10. P. 668-676. DOI: 10.1061/(ASCE)IR.1943-4774.0000337.

Статья получена 20.11.2016 г.

Churuksaeva V.V., Starchenko A.V. NUMERICAL INVESTIGATION OF A TWO-PHASE FLOW OF FLUID WITH LIGHT PARTICLES IN OPEN CHANNELS Tomsk State University Journal of Mathematics and Mechanics. 6(44). pp. 88-103

DOI 10.17223/19988621/44/8

A mathematical model and a computational method for a numerical investigation of the two-phase turbulent flow in an open channel are performed. The solid particles with a density close to that of water were considered as a dispersed phase. The model is based on the flow depth-averaged equations of mechanics of interacting and interpenetrating continuums in a hydrostatic approach. A turbulent closure of the model is implemented with the application of the k — e turbulence model modified by Pourahmadi and Humphrey (1983) to consider the influence of particles on the turbulent structure of the flow.

The numerical method proposed for solving equations of the model is based on the elimination algorithm and explicit-implicit time approximation.

An unsteady turbulent flow in a 180-degree bend flume with polypropylene particles modeling the ice was computed and the results were compared with those of Urroz and Ettema (1992). It was found that the mathematical model and the computational method proposed predict accurately both the velocity field and distribution of the particles in the channel.

The influence of the dynamic parameters of dispersed phase on the turbulent structure of the flow was investigated by conducting the calculations of the flow in an open channel with a 90-degree bend. It was revealed that the structure of a two-phase flow is most affected by the size and shape of the particles.

Keywords: mathematical modeling, two-phase flow, double-speed continuum, shallow water approximation, k — e turbulence model, ice particles, finite volume method.

CHURUKSAEVA Vladislava Vasilievna (Tomsk State University, Tomsk, Russian Federation) E-mail: chu.vv@mail.ru

STARCHENKO Alexander Vasilievich (Doctor of Physics and Mathematics, Tomsk State University, Tomsk, Russian Federation) E-mail: starch@math.tsu.ru

REFERENCES

1. Nigmatulin R.I. (1991) Dynamics of multiphase media. New York: Hemisphere Publishing Corporation. 375 p.

2. Bubenchikov A.M., Starchenko A.V. (1998) Chislennye modeli dinamiki i goreniya aerodis-persnykh smesey v kanalakh [Numerical models of the dynamics and combustion of aerodis-perse mixtures in channels]. Tomsk: Tomsk University Press. 236 p.

3. Churuksaeva V.V., Starchenko A.V. (2015) Matematicheskaya model' i chislennyy metod dlya rascheta turbulentnogo techeniya v rusle reki [A mathematical model and numerical method for computation of a turbulent river stream]. Vestnik Tomskogo Gosudarstvennogo Universiteta. Matematika i mekhanika - Tomsk State UniversityJournal of Mathematics and Mechanics. 6(38). pp. 100-114. DOI 10.17223/19988621/38/12.

4. Rodi W. (1980) Turbulence models for environmental problems. In: Prediction Methods for Turbulent Flow, von Karman Inst. New York: McGraw-Hill. pp. 276-281.

5. Launder B.E., Spalding D.B. (1974) The numerical computation of turbulent flows. Computer Methods in Applied Mechanics and Engineering. 2(3). pp. 269-289. DOI 10.1016/0045-7825(74)90029-2.

6. Gidaspov D. (1994) Multiphase flow and fluidization: continuum and kinetic theory descriptions. Boston: Academic Press. 457 p.

7. Patankar S. (1980) Numerical heat transfer and fluid flow. CRC Press, 214 p.

8. van Leer B. (1979) Towards the ultimate conservative difference scheme. V-A second order sequel to godunov's method. Journal of Computational Physics. 32. pp. 101-136. DOI 10.1016/0021-9991(79)90145-1.

9. Cada M., Torrilhon M. (2009) Compact third-order limiter functions for finite volume methods. Journal of Computational Physics. 228. pp. 4118-4145. DOI 10.1016/j.jcp.2009.02.020.

10. Cea L., Puertas J., Vazquez-Cendon M.E. (2007) Depth averaged modelling of turbulent shallow water flow with wet-dry fronts. Archives of Computational Methods in Engineering. 14(3). pp. 303-341. DOI 10.1007/s11831-007-9009-3.

11. Hou J., Simons F., Mahgoub M., Hinkelmann R. (2013) A robust well-balanced model on unstructured grids for shallow water flows with wetting and drying over complex topography. Computer Methods in Applied Mechanics and Engineering. 257. pp. 126-149. DOI 10.1016/j.cma.2013.01.015.

12. Urroz G.E., Ettema R. (1992) Bend ice jams: laboratory observations. Canadian Journal of Civil Engineering. 19. pp. 855-864. DOI 10.1139/l92-097.

13. Tsai W.F., Ettema R. (1994) Ice cover influence on transverse bed slopes in a curved alluvial channel. Journal of Hydraulic Research. 32(4). pp. 561-581. DOI 10.1080/00221686. 1994.9728355.

14. Han S.S., Ramamurthy A.S., Biron P.M. (2011) Characteristics of flow around open channel 90° bends with vanes. Journal of Irrigation and Drainage Engineering. 137(10). pp. 668676. DOI 10.1061/(ASCE)IR.1943-4774.0000337.