АЛГОРИТМЫ ПРЯМОГО ЧИСЛЕННОГО МОДЕЛИРОВАНИЯ ДИНАМИКИ ДИСПЕРСНОЙ ФАЗЫ ПРИ ОБТЕКАНИИ ТЕЛА ЗАПЫЛЕННЫМ ПОТОКОМ Текст научной статьи по специальности «Физика»

УДК 532.529.5

Алгоритмы прямого численного моделирования динамики дисперсной фазы

при обтекании тела запыленным потоком1

Д.Л. Ревизников, А.В. Способин

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

Введение

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

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

1 Работа выполнена при поддержке РФФИ (проект № 05-08-01478-а)

• с высокой точностью определять характеристики дисперсной фазы в заданной точке в конкретный момент времени;

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

• исследовать переходные процессы, возникающие при входе объекта в запыленное облако;

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

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

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

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

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

й V.

т

р _

= F

I

р йг

й® р -_р_т

р йг _ и

где

I.

- масса, момент инерции, скорость и угловая скорость частицы,

Р Р внешние силы,

0),

- момент внешних сил. В качестве внешних сил, приложенных к частице,

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

и сила Магнуса

вызванная

вращением частицы. На угловую скорость оказывает воздействие вращающий момент

Р

D

М

о

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

определяется выражением

2

где - радиус частицы,

ГР

- плотность и скорость газа. Коэффициент сопротивления

Р V

иа 9

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

^ (Re рМр)

настоящей работе соотношением Хендерсона [6].

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

РМ = П \С-Рд

1 VxV -со х(Т -V )

И д р \ д Р1

Коэффициент определяется следующим образом [7], [8]:

с =с (Ие Де )

— — \ р' — )

Сс =

Ие

-< 0.45

0.45

Ие р

Ие,

■ + |1-0.45—р -ехр(-0.05684 Яе?; Ке0п'3), —->0.45 Иес - р )' Иер

\ - р

Выражение для вращающего момента имеет вид

Т = Грс,р (1 vxv -v )|1 vxv -v |

- 2 9\2 д р 2 9 р

Выражение для коэффициента момента

с=С1 (Ке-)

заимствовано из работы [9]:

•< («е- )=

64 п Ие.,

Ие <32

Здесь

М.

12.9 . 128.4 „ „

—-, Ие ?>32

Ие0 5 Ие? ' -

со ш

- число Маха, ,

Неп Кем

р о

- числа Рейнольдса для относительного поступательного

и вращательного движения соответственно:

1

М 4^! Ке_2грр.1'ур4Ухт9-5р|

Р ^ Р Ря Ке„_ 2

Ря

где - газовая постоянная, - показатель адиабаты, - коэффициент вязкости газа,

" ' Ря

определяемый по формуле Саттерленда, - температура газа.

Т

Для определения параметров частицы после соударения с другой частицей используется модель твердых сфер [3]. В ее основе лежат уравнения для импульса и момента импульса системы двух частиц:

т1(Т1-Т1°1)_Т ,

т2( Т 2—Т(20) 1_—п

(Т 2—Т20))_-

1 (о) 1 — ))10))_ г 1 пх7

(ы 2—Т 20))_г 2п х7

12 (Т2 — 4 )

где ^ - единичный вектор, направленный из центра масс частицы 1 в центр масс частицы 2,

п г1

и - радиусы, и - массы, и - моменты инерции первой и второй частиц

г2 т1 т2 11 12

соответственно. Индексом обозначены параметры частиц до столкновения. Данная система

(0)

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

-г, Т Хп + Гп

С'0)_С(0)+Г 1 т10)хТ+г 2 Т'0)хп

Здесь , - относительные скорости центров масс частиц до и после

£(°)_п 10)—Т0) 5_п1—п2

удара. Тангенциальная компонента относительной скорости ¿(о)_£(„)—(¿(о-п) П_£(о) — (£(о-п) п+ Гх "10хп+ г2 т'0)хп '

Импульс может быть представлен в виде суммы нормальной и тангенциальной составляющих

где

V = J п + Jtt

п I

¿(0)

. Коэффициенты восстановления и трения

е г

вводятся как

п-о=-в(т-с(0Г ^=fJn

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

множества решений. Если выполнено соотношение

п-£(01 2 1

то скольжение

й)! 7 f (1+ е)

продолжается в течение всего процесса взаимодействия, и параметры частиц после удара равны:

V , = у Г'-(V - V) (т-б(о)) (1 +е)

1= у 1

тп

т1+т2

V2=V20)+(п-fV)(v•G(o')(1+ е) —1

т1+т2

2

(V•G(0')(пхV)[(1+е)- т 1 1 2 г / т1+т2

1

V, ^'--(v•G(0))(VXV)Г (1+е)- т 2 2 2 г2 т1+ т2

В противном случае скольжение прекращается, и относительная скорость точки контакта

становится равной нулю, поэтому

t ^ 7} т1+т2 а1

, и решение:

►,(0).

у 1= у 1 -

►,(0)

у 2= у 2 +

(1+е) (v•G(o)) п+7!GС0)|

(1+е)(V•G(0') п+2Й'|

тп

т1+т2

т,

т1+т2

Тl=тi0i-тг- ¡¿с; '| (п х V)

тп

1

т1+т2

1

- -(0) 5 ,яф),/^-\ т а = ы 2'--Юу |( пх г )-

2 2 7 г 2' сН1 т1+т2

При расчете взаимодействия частицы с поверхностью обтекаемого тела используется модель твердых сфер [3]. В ее основе также лежат уравнения для импульса и момента импульса. В процессе взаимодействия выделяют период сжатия и период восстановления. Скольжение частицы

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

- касательная к поверхности в точке удара, ось ортогональна плоскости .

х—г у х—г

Если скорость частицы удовлетворяет условию , то решение выглядит

(0)

< 2 |Т(0)| (е+1)

следующим образом:

(0) '

Г

Ыу _ы(у°) _ ^ —--

г г

В противном случае, если верно , то параметры частицы после отражения

(0)

2 V

—ттте+гу <0

от поверхности определяются формулами: vх_vХ0)+ех( (е+1) VУ) ,

vу _—evу

У) '

VZ_ V ? ^ (е +1) V У0)

)х _о(х0) — ^ (е+1) $>

-„( 0) '

Ыу _Ыу

)_4°'+(е+1) vУ01

Здесь , - направляющие косинусы проекции скорости до удара в плоскости

£х £г

Коэффициент восстановления не является постоянным и определяется соотношением [10]

е ( у ) =

Здесь:

1-( 1 -е 0)

/ \

I \1 у_

/0,

-1/

5 > у < у0

^ у> у0

\ и/

модуль нормальной составляющей скорости в момент соударения с поверхностью, а

и определяются характеристиками материала.

е0 у0

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

уу е^У

(0)

у

X'

4°

(0),

ууХ'ет+— " г (ет-1).

у (х0) ет - 2 — (0) г,

XI 7 '

— =

5 у(х°)(ет-1) + 5 —(0)

'X (ет_

уу е

2 г 2

(0) е 2 ^ + 2 —( 0)

г 7

ет- 5

в< в в > в

в< в1 в > в

Здесь ось направлена по нормали к поверхности в точке удара, ей ортогональна, ( )

У X у10)

.,(0)

со

(0)

-компоненты скорости и угловой скорости частицы до и после

соударения,

коэффициенты восстановления нормальной и тангенциальной компонент

ем ет

импульса, - радиус частицы, - угол между

г в у

,, и осью

(0) X

в'

- критическое

значение угла.

у

е

0

у

3

x

У

Особенности реализации модели

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

систему координат . Ось совпадает с осью цилиндра, ось направлена вдоль

Охуг Ог Ох

течения, ось им ортогональна. Определим число частиц в двумерной области

Оу

, соответствующее объемной концентрации примеси . Рассмотрим объем

(ха'хБ)х(УА'УБ) Су

, равномерно заполненный примесью. Объем, занимаемый примесью,

( ха'хБ)х( уА'уБ)х( 0 М

равен . Вследствие равномерного распределения примеси все сечения

ур_су ( хб—ха )( уб — уа )н

данного объема плоскостями, параллельными равноправны. В каждом сечении примесь в

Оху

среднем занимает площадь . Поскольку частицы распределены в

5р_-у _СУ ( хб — ха)( уб — уа ) объеме равномерно, их центры масс в общем случае не лежат в плоскости сечения. Средняя

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

1 г 2

5 _— I п (г2 — г2) йг _ —пг2

2 г — 1 1 3

—г

число частиц радиуса в области при объемной концентрации примеси

г ( ха'хб)х( уа'уб)

равно . В полидисперсной примеси каждому сорту частиц

Су 5р _ ЗСу!ХБ—УА)

5 2 пг 2

соответствует объемная концентрация и число частиц в

Сук _ЯкСу К

I Як_1

к _1

рассматриваемой области .

„ _3 ЯкСу(хб—ха)(уб — уа)

^ к _ г-, 2

2 пг

Расчет производится с шагом по времени . Для решения уравнений движения применяется

т

метод Рунге-Кутты пятого порядка точности. После вычисления перемещения частиц за интервал

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

(гк—1 'Л ]

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

хк )_акх' 2 + Ькх* + Скх

Ук(I )_аку^2+ Ьку' +Ску

Условием соударения частиц и служит равенство расстояния между их центрами масс

1 j

сумме радиусов . Представленное уравнение четвертой

(х; (Г) — х, (Г ))2+( у; (Г) — yJ (' ))2_(г. + г] )2

степени решается путем сведения его к уравнению третьей степени с помощью специального вида замены переменной [12]. Для каждой пары находящихся в зоне досягаемости друг для друга за

интервал частиц проверяется возможность столкновения. Формируется единая очередь

('к—1 ]

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

частицей , влечет за собой расчет ее параметров после удара, аннулирование всех 1

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

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

Результаты

В данном разделе иллюстрируются возможности разработанного программно-алгоритмического аппарата.

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

табличных данных [13,14]. Материал частиц - песок с плотностью 2500 , диаметр частиц

кг/м

примеси 5, 10 и 50 мкм. Все частицы в каждом из экспериментов считались однородными шарами одинакового размера и плотности. Объемная концентрация дисперсной фазы в области

невозмущенного течения составляла и .

10-5 10-4

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

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

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

построении графиков характеристик применялось осреднение по пространству с шагом ,

1°

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

осреднения в зависимости от режима варьировался от до с. Учитывались удары со

210-4 10-2

скоростями, превышающими 1 м/с.

На рис.1 представлены распределения частиц диаметром 5 мкм в различные моменты времени сразу после попадания тела в пылевое облако. Результаты получены по бесстолкновительной

модели при концентрации примеси в области невозмущенного течения . Видно постепенное

10-5

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

Рис. 1. Эволюция распределения частиц диаметром 5 мкм без учета соударений и вращения.

На рис.2 представлены установившиеся распределения частиц, полученные по трем различным моделям. Учет взаимодействия частиц друг с другом приводит к размыванию границы рассматриваемой и повышению концентрации примеси вблизи поверхности. Вращение частиц способствует некоторому расширению зоны повышенной концентрации, однако существенного вклада в картину распределения не вносит.

Рис. 2. Установившееся распределение примеси частиц диаметром 5 мкм в режимах: а) без соударений и вращения; б) с учетом соударений; в) с учетом соударений и вращения.

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

Из графика на рис.3 видно, что учет соударений частиц приводит примерно к двукратному

росту интенсивности ударов о поверхность тела уже при объемной концентрации примеси .

10-5

Повышение концентрации до влечет многократное увеличение числа ударов частиц о тело.

10-4

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

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

Учет вращения частиц вызывает дальнейшее ослабление теплового воздействия практически на всей поверхности цилиндра.

1,8-1 1,6 1,4 1,2 1,00,80,6 0,4 0,2 0,0

■ ■ без учета соударений

----с учетом соударений

-о—о- с учетом вращения

C =10

V0

D = 5 мкм

2018161412 ■ 108 6 4 2 0

■ ■ без учета соударений Су0-10

с учетом соударений р = 5 мкм -о—о- с учетом вращения

1—1—г -80 -60 -40 -20

а, град

"Г 0

а, град

г^ I—1—г 40 60 80

Рис. 3. Интенсивность ударов о тело частиц диаметром 5 мкм.

0,30-, 0,250,20->§ 0,15-

Z >

0,10 0,05

0,00

■ ■ без учета соударений Су0= 10

с учетом соударений р - 5 мкм

с учетом вращения

-i—I—|—I—|—I—|—I—|—I—|—I—|—I—|—I—|—

■80 -60 -40 -20 0 20 40 60 80

а, град

0,30 0,25 0,20 >Z 0,15 ■

Z >

0,10 0,05 0,00

■ ■ без учета соударений

--- с учетом соударений

-о-о- с учетом вращения

C =10-4

V0

D = 5 мкм

-i—I—|—I—|—I—|—I—|—I—|—I—|—I—|—I—|—

-80 -60 -40 -20 0 20 40 60 80

а, град

Рис. 4. Средняя нормальная составляющая скорости частиц диаметром 5 мкм в момент соударения с поверхностью.

0,06-1

0,05-

0,04

g° 0,03 О

0,02 -

0,01

0,00

■ ■ без учета соударений

----с учетом соударений

с учетом вращения

C =10

V0

D = 5 мкм

И 1 I 1 I 1 I 1 I 1 I 1 г -80 -60 -40 -20 0 20 40 60 80

а, град

0,06-,

0,05-

0,04

0° 0,03

°

0,02 -

■ ■ без учета соударений Суо-10

с учетом соударений р = 5 мкм н>ч> с учетом вращения

0,01

°,°0 -I-1-1-Г"1-1-1-1-1-1-1-1-1-1--1-1-1

-80 -60 -40 -20 0 20 40 60 80

a, град

Рис. 5. Удельная мощность воздействия на тело примеси частиц диаметром 5 мкм.

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

1. Нигматулин Р.И. Основы механики гетерогенных сред. М.: Наука, 1978, 336 с.

2. Стернин Л.Е., Маслов Б.Н., Шрайбер А.А., Подвысоцкий А.М. Двухфазные моно и полидисперсные течения газа с частицами. М.: Машиностроение, 1980, 172 с.

3. С.Т. Crowe, M. Sommerfeld, Y. Tsuji. Multiphase flows with droplets and particles. CRC Press LLC, 1998, 471 р.

4. Гилинский М.М., Стасенко А.Л. Сверхзвуковые газодисперсные струи. М.: Машиностроение, 1990, 176 с.

5. Tsirkunov Yu. M. Gas-particle flows around bodies - key problems, modeling and numerical analysis. // Proc. Fourth International Conference on Multiphase Flow (Ed.: E. Michaelides), May 27 - June1, 2001, New Orleans, LA, USA. - CD ROM Proc. ICMF'2001, paper No. 609, 31 p.

6. Henderson C.B. Drag coefficients of spheres in continuum and. rarefied flows. // AIAA Journal Vol. 14, No. 6, June 1976. P. 707-708.

7. Rubinow S.I., Keller J.B. The transverse force on a spinning sphere moving in viscous fluid. // J. Fluid Mech. 1961. V. 11. Pt. 3. P. 447-459.

8. Oesterle B., Bui Dinh T. Experiments on the lift of a spinning sphere in a range of intermediate Reynolds numbers. // Experiments in Fluids. Vol. 25, No. 1, June 1998. P. 16-22.

9. Dennis S.C.R., Singh S.N., Ingham D.B. The steady flow due to a rotating sphere at low and moderate Reynolds numbers. // J. Fluid Mech. 1981. V. 101. Pt. 2. P. 257-280.

10. McNamara S., Falcon E. Simulations of vibrated granular medium with impact-velocity-dependent restitution coefficient. // Physical Review E 71 2005. 031302-1 - 031302-6.

11. Циркунов Ю.М., Панфилов С.В., Клычников М.Б. Полуэмпирическая модель ударного взаимодействия дисперсной частицы примеси с поверхностью, обтекаемой потоком газовзвеси. // Инж.-физ. журн. 1994. Т. 67. № 5, 6. С. 379-386.

12. Подвысоцкий В. Общий аналитический метод решения алгебраических уравнений четвертой степени. // Электронная библиотека «Наука и техника». - http://www.n-t.org/tp/ns/oam.htm (17.11.2006)

13. Любимов А.Н., Русанов В.В. Течения газа около тупых тел: в 2 ч. Ч. 1. М.: Наука, 1970. - 287 с.

14. Любимов А.Н., Русанов В.В. Течения газа около тупых тел: в 2 ч. Ч. 2. М.: Наука, 1970. - 379 с.

Сведения об авторах

Ревизников Дмитрий Леонидович, профессор кафедры вычислительной математики и программирования Московского авиационного института (государственного технического университета), д.ф.-м.н; e-mail: reviz@k806.mainet.msk.su; контактный телефон: 158-48-94;

Способин Андрей Витальевич, аспирант кафедры вычислительной математики и

программирования Московского авиационного института (государственного технического

университета);

e-mail: spise@inbox.ru;

контактный телефон: 777-00-55.