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

УДК 004.021:536.46

DOI: 10.18698/1812-3368-2018-6-75-91

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

Т.Н. Романова Д.А. Ягодников Г.А. Щетинин

rtn@bmstu.ru daj@bmstu.ru gashetinin@yandex.ru

МГТУ им. Н.Э. Баумана, Москва, Российская Федерация

Аннотация

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

Ключевые слова

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

Поступила в редакцию 09.04.2018 © МГТУ им. Н.Э. Баумана, 2018

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

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

Моделирование горения ПМ, например алюминия, основывается на фундаментальных законах переноса энергии, количества движения и массы реагентов с учетом наличия объемных источников тепловыделения и массы [8, 9]. Однако сложности объекта моделирования, обусловленные наличием газовой и дисперсной фаз, протекающих как в кинетическом, так и в диффузионном режимах химическими реакциями, а также распределение частиц ПМ по размеру приводили к необходимости принятия нескольких допущений. Например, допущения о равенстве скоростей газовой и дисперсной фаз [10], а также о монофракционном составе газодисперсной реагирующей системы [11]. Кроме того, большинство математических моделей реализовывалось в одномерной (по пространственной координате) постановке ввиду сложности вычислений [12].

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

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

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

- отражение реального полидисперсного распределения частиц;

- учет различия скоростей и температуры газа и частиц;

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

Для отражения в модели перечисленных требований можно использовать перспективный метод, связанный с применением аппарата статистической физики и, в частности, функции плотности распределения вероятности (ПРВ) одной из нескольких физических величин [13-15]. Основное достоинство указанного метода состоит в возможности учета реального полидисперсного распределения частиц ПМ, всех режимов горения (ламинарного и турбулентного), а также влияния турбулентных пульсаций различных параметров газовой взвеси на осредненные значения скоростей химических реакций и тепловыделения. Такая постановка задачи соответствует сложившимся представлениям о взаимосвязи термодинамических, кинетических и газодинамических процессов.

В общем случае аргументами функции ПРВ, кроме пространственных координат и времени, являются термодинамические параметры газовой и конденсированной фаз. В настоящей работе для описания процесса горения в движущейся аэровзвеси частиц алюминия используется функция совместной плотности распределения, аргументами которой приняты время I, декартовы координаты х, у, радиус Гк и температура Тк частицы алюминия: Р(£ х, у, Тк, п). Будем пола-

гать, что Р(£ х, у Тк, Гк)йТкйтк определяет вероятность того, что для частицы, находящейся в момент времени I в координатах (х, у), значения гк, Тк принадлежат к интервалам гк е (гк, гк + йтк) и Тк е (Тк, Тк + йТк) соответственно. При этом, используя понятие и уравнение баланса функции ПРВ, контролируемое условием нормировки, необходимо ввести дифференциальное уравнение баланса счетной концентрации частиц Пк:

дПк дПкР дПкР $2ПкР

—— +Я—— йткйТк +|Г—— йткйТк - Л--—— йткйТк -

8х п дх п ду п 2 дх2

,, D д2пкУ г, D г, D 52nkP

-Я 1 lXydrkdTk 7W k k "ЯT

dnkP dnkP dnkP D (d2nkP 82nkP 82nkP

dz

8x

8y 2 ^ 8x2 8x8y 8y8x dan P , 8fn P = 0

drkdTk = 0;

sn?"

8y2

(1)

-vr n - (2)

8Tk 8rk

где nk — счетная концентрация частиц; D — коэффициент диффузии частиц ПМ; &k — скорость изменения температуры частицы, K/с; fk — скорость изменения радиуса частицы, м/с.

Граничные условия для уравнений (1), (2) зададим следующим образом:

8щ

x = 0:

= 3WMДвозТ ;

4nrimaxmo2 ' P = Y(rk)y(Tt -Tkо);

x — x«

8x SP

- 8x

■ = 0; = 0;

y = 0:

8nk 8y

= 0;

— - 0;

8y

y=yn

^L = 0;

8y

d-P = 0.

8y

Здесь Wai, mO стоянная для воздуха; y(rk) =

— относительные массовые концентрации; Евоз — газовая по-

(rk -0,5)2

— начальная функция распреде-

1

0,18

0,зТ2л

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

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

др

-div(up) = -£ >Gj;

dz '

p = (1 - z )pRT;

,5T

pC--h upC grad T = div((A, + X т) grad T) + Q;

dz

dmi

+ pu grad mi = div(p(D + DT) grad mi) ± Gi,

(3)

(4)

(5)

где р — плотность газа; т — время; u — скорость газа; p — давление газа; 2 — относительная массовая концентрация компонентов к-фазы; Я — универсальная газовая постоянная; Т — температура газа; С — теплоемкость; X, Хт — коэффициенты молекулярной и турбулентной теплопроводности газа; О — источниковый член, определяющий поток теплоты в газовую фазу; Ш{ — относительная массовая концентрация компонентов; Б, — коэффициенты молекулярной и турбулентной диффузии газа; 0} — источниковый член, определяющий массовую скорость образования или исчезновения к-го компонента газа.

Граничные условия для уравнений (3)-(5) зададим следующим образом:

x = 0: <

x — Ха

^ = 0;

5x

f=* dx

U — UH ;

T = T ;

н

Wo2 = 0,23; wa = 0,77;

^ = 0; öy

y = 0:

dml dx

= 0;

f=0;

oy

У = Ут

5ml öy

= 0;

^ = 0; öy

f=0 Oy

dml

dy

= 0.

Здесь ин, Тн — начальные значения компонентов скорости и температуры газа; модуль скорости газа при у = 0 и у = утж изменяется согласно условию скольжения. Зададим плотность р и теплоемкость С газа, коэффициенты молекулярной теплопроводности Л. и молекулярной диффузии газа Б как функции температуры газа [15]:

р = _Рат^; я = 26,2.Ю-3 +(Т-300)-7,6-10~5; С = 10~3 + 0,87(Т-300);

ЯвозТ

D = 0,178 -10~4 -

4T 273

где ратм — атмосферное давление.

Источниковые члены 0} и О в уравнениях (3)-(5) определяются с учетом конвективного и радиационного механизмов передачи теплоты и в предположении протекания одностадийной химической реакции

А1 + (3/4)02 ^ (1/2)АЬ0э + АН с помощью функции ПРВ следующим образом:

4

О = АН [Г 4то-к2ркПк/кРйТкйтк - |Г-то-к2рСкРк&кЩ/кРйТкйтк;

3 (6)

0 = Кт0 Ц 4лгк2ркПк /к РйТкйтк.

Здесь Кт0 — массовый стехиометрический коэффициент.

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

( ¿т- а ■ 2п п_1

/к = —Т =-01,^02 пТкп 1, (7)

йх Пок р (Т) )0,2

где пок — концентрация эффективного окислителя; п — показатель степени; Т0 — начальная температура в камере сгорания.

Изменение температуры частицы за единицу времени рассчитывается по следующей формуле [15]:

йТк _О_

®к =—- =---. (8)

й Ащгг2 I ТкС2р2

4%т- у~2 +51С1р1

Здесь индекс «1» обозначает параметры, относящиеся к продукту химической реакции (окисление), а индекс «2» — к исходным компонентам.

Начальные значения параметров в камере сгорания при т = 0 (обозначим их индексом «0»): давление р0 = 0,1 МПа, температура Т0 = 300 К, скорость движения газа Ы0 = 1 м/с, относительные массовые концентрации т02 = 0,23, тА1 = 0,77. В качестве горючего рассматривается порошкообразный алюминий. Воспламенение частиц осуществляется горячим газом, который находится во всем объеме Х0< х <Хж и имеет постоянную температуру Тв = 2300 К (частицы достигают такой температуры на фронте пламени). Скорость распространения фронта пламени рассчитывается по массовой скорости выгорания компонента, находящегося в недостатке [15]:

1

Wf =

Р0 J Go2 dX' 0,89а

а < 1; а>1,

Р0 J GkdX

где а — начальное значение коэффициента избытка окислителя.

Геометрическое описание расчетной области. Моделируются процессы воспламенения и горения аэровзвеси частиц алюминия в камере сгорания с прямоугольным сечением. Расчетные области процесса приведены на рис. 1.

Г\

Фронт пламени

Уо о'ШШШМШШШШв. х х0 = о ж

Рис. 1. Геометрия расчетной области и общая схема процесса

Поставленная задача решается численно в двумерной постановке с декартовыми координатами, где x0 = 0 — координата на левой границе, на которой осуществляется подача частиц в рабочую зону (камеру сгорания). Частицы при подаче обладают характеристиками, начальные значения которых поддерживаются постоянными: ukн (укн, ) — начальная скорость частицы; щн , wli) — начальная скорость газа; Тк — температура частицы; гк — радиус частицы.

На верхней (у = утах ) и нижней (у = 0) стенках применяются условия скольжения: проекция скорости (как газа, так и частицы) на ось ОХ становится равной по модулю скорости на предыдущем слое, а проекция скорости на ось ОY — нулю, теплообмена со стенкой и внешней средой не происходит (адиабатическая стенка), относительная массовая концентрация не меняется, счетная концентрация частиц, как и значение функции ПРВ, не изменяется. На правой границе (х = x<ю) применяются «мягкие» граничные условия.

В начальный момент времени частицы находятся на левой границе x0 = 0, распределены равномерно вдоль вертикальной координаты Y и обладают заданными начальными значениями температуры Tk, радиуса Гк и скорости щ. При нагреве горячим газом, находящимся во всем объеме камеры сгорания, и за счет тепловыделения в результате химической реакции параметры частиц меняются, согласно уравнениям (7), (8).

Численный метод решения. Дифференциальные уравнения (1)-(5) аппроксимируем разностными уравнениями первого порядка точности по времени и первого порядка точности по пространству. Для этого введем равномерную разностную сетку по осям координат ОХ ( = 0,...,Ы), OY (] = 0,...,М), по временной оси От и по осям параметров частиц — Оп (д = 0,...,дтах) и ОТк ($ = 0,..., 5тах). После тестирования определено, что наиболее качественные результаты получаются при значениях N = 24, М = 24, дтах = 14, 5тах = 14.

Уравнение неразрывности

ф f dv dw Л ф ф

+ Р

дг

дх ду ^ дх ду

где V — проекция скорости газа на ось ОХ; V — проекция скорости газа на ось ОУ. Далее переменные, которые взяты с предыдущего временного слоя, обозначены волнистой чертой. В результате аппроксимации левой разностной производной получаем линейные уравнения по каждой координате скорости

_1,) + BxVit) + Сх = 0; ¡-1 + Ву) + Су = 0,

где

2Pi,j -Pi-i,j C _Pi,j - Pi,j , i v G

7 ' - + ^ jGj ;

2x 2 J

hx

Ax = ^; = hx

A _~Pi,j ; B _ 2Pi,j -Pi,j-i ; C _Pi,j - Pi,j , 1 v G

— : ; Dy — - ; Cy — - . *

W

h

y

h

y

y

2T

jj

— прогоночные коэффициенты. Уравнение баланса теплоты

рС(Т)

(дТ дТ дТ') д Л дТ , . — + v — + w— = — I X— 1+—

V Q.

дг dx dy) dx \dx) dy ^ ду В результате аппроксимации получаем разностные уравнения

Т

¿г. *

-рi,jC(T,j ) +Т,j pi,jCTj)

V

1+Vij+Wj

T hx hy J

(

"Xi+i,j j ^i,j+i j

hx

h

2

y у

Л

h

hy 2

AtT-1, j + DT j_i + CtT, j + DT j+i + ETT+1, j = - Ft . Здесь прогоночные коэффициенты имеют следующий вид:

+Т-1,j | ~Pi,CTj )Vi,j -— %itj —■ ] + T,j-i "Pi,jC(Ti,j )wi,j -— %itj —

hx hx J у hy hy

T T

л 1l+i,j 1 1t,j+i n _ n.

i+i, j—^—^i, j— Q - 0;

(9)

Pi, jC(Ti, j )vi, j Xi> j pi, jC(Ti, j )wi, j Xi> j at -------—-; dT -----■

hx

ct =p i, jC(T j)

hx

i+VViii I w ,j

x hx h

hy

hy'

^i+i,j + ,j , ^i,j+i + ,j

y У

hy2

D = ^,j+i . £ = ^i+i,j . ^ = Pi,jC(Ti-,jj _ Q;

hy hx T

^i,j -

,j +^i-1,j ; y j ,j-1; у _ ^i,j +^i+1,j; T - ^i,j ,j+1 ; ^i,j _ ~ ; ^i+1,j _ ~ ; ^i,j+1 _

2 2 2 " 2 Для решения системы линейных уравнений (9) использован метод Гаусса для пятидиагональной матрицы, в которой записаны прогоночные коэффициенты. Уравнение сохранения массы

(дщ дт, дт, Л д ( дщ ^ д ( дщ

дг dx dy ) dx v dx ) dy ^ dy заменяем разностным уравнением

m

',)

m

',)

где

^ ( 1 V,j Wi,j ^ -

Pi,j -+-TL+T2

pi +1, j Di+1 ,j +Pi,jDi,j . Pi,j+1Di,j+1 +Pi,jDi,j

hx

- mi_1, j | -LL L m, j_11 - ml

hhx ^hv

hy2

Л

hy hy J

Wi+1, jPi+1, jDi+1, j mi, j+1Pi, j+1Д, j+1 _ ----+Gi - 0;

hx 2

hy2

Ammi—1, j ^ Bmmi, j—1 ^ Cmmi, j ^ Dmmi, j+1 ^ Emmi+1, j Fm,

A _ Pi,jyi,j Pi,jDi,j; B __Pi,jW,j Pi,jDi,j

hx ^hx

hy hy

Cm ~ Pi,j

1 + ViL +WlL

у T hx hy J

Л -К

Pi+1, jDi+1, j +Pi, jDi, j Pi, j+Di, j+1 +Pi, jDi, j

hx

hy2

Dm —

Pi, j+Di, j+1

hy2

; E _ pi+1,jA+1,j; F +G |;

> ^m ~ 1 ^ > rm ~ | F — ^i | >

(10)

_ Pi, j +Pi-1, j _ Pi ,j +Pi, j-1 _ Pi ,j +Pi+1, j _ Pi, j +Pi, j+1 Pi,j ~-~-; Pi,j ~-~-; Pi+1,j ~-~-; Pi,j+1 _-~-

Di, j =■

-; Di+1,j =■

-; Di,j+1 =-

- В,,} + В,_1,; — _ Вц+р, ^_1 - _ В,,) + В,+1,1 - _ Д,) + В,^+1

-; в,, 1 —-

2 2

— прогоночные коэффициенты.

Для решения системы линейных уравнений (10) использован метод Гаусса для пятидиагональной матрицы, в которой записаны прогоночные коэффициенты. Уравнения, описывающие источниковые члены:

_ кТ кг „ кТ кг

Тк + Т гк+Т Тк + +кТ 4

О = АИ X Е 4лгк2ркИк/кР- Е Е - шкРСкРк®кПк/кР;

кТ кг _ кТ кг 3

Тк —— гк-— Тк —— гк- —

2 2 2 2

„ кт hr

Tk+hT rk+IT

G = Km0 E E 4^rk2pknkfkP.

Tkrk 2 2

Уравнение баланса счетной концентрации частиц в дисперсной фазе:

-йткйТк -

дПк dkT + tf^pdrkdTk-Ц D ß2nkР

dz

dy

2 dx2

r, D c^nP rr D <Ы?ЙЙТ rr D S2nkP

2 dxdy JJ 2 dydx JJ 2 Sy2

Введем верхние индексы «q» и «s» по координатам rk и Tk:

drkdTk = 0.

nk ;i,j nk;i,j

rmax Tmax

-z z

ro To

( (

nk ;i, j

V V

2pq, _2pqj -pq^ d(ps,+2pq,js -pfj5!.j)

2hx

d (4pq:; - p^, j -

^) d (pq,^+2pq,js - p;)

2hxhy

2hy2

(

nk ;i-1, j

-pq,j -Dpq,js d(2pq:si,j -pqj -pj

2hxhy

+ nk ;i ,j-1

h,x 2hx

-pq,j -Dpq,js d(ipq:;_i-pq^H-pq,;)

л

hv

2hy2

2hxhy

-d (pq:h} - pqj) -d (pji - pqj)

nk ;i+1, j-—2-+ nk ;i, j +1 —

2hx 2

-D

+ nk ;i-1, j-1 "

/pq,s _ pq,s \ ^pi, j-1 pi-1, j)

2hx hy

2h2

Л

hrhT = 0;

Annk ;i ,j + Bnnk ;i-1, j + Cnnk ;i, j-1 + Dnnk ;i+1,j + Ennk ;i ,j +1 + Fnnk;t -1, j-1 Gn.

Здесь прогоночные коэффициенты 1

An — '

rmax Tmax z z ro To ^ opq,s _ z'ri, j P«q-1, j 2

hx

d (щ; _ pq,s i-1, j _ pq,s \ pi, j-1)

2hxhy

rmax Tmax = Z Z ro To ' pq,s ^ j - DPq,s

hx V i ' 2hx 2

rmax Tmax = 5: z ro To ' pq,s ^ j - DPq,s l ', j

hy V ' 2hy2

hy

2hx 2

2hy2

hrhT;

2hx hy

hhT;

2hx hy

hrhT;

(11)

rmax Tmax

Dn = Z Z

-0 To

-D ( p*î, , - pq; )

2h*

rmax Tmax

En = z z

П5 To

-D

(pq,s pq, s\ ^ ¿¿+1 P ¡,j j

2hy2

hh ;

hh ;

rmax Tmax

Fn = z z

ro To

-D

(pq,s _ pq,s \ ¡, j-i p ¡-1, j )

2hx hy

л

h-hT ; Gn =■

-nti, j

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

ônkp + ônkp + ônkp D Г ô2nkp + d2nkp + d2nkp + ô2nkp дг dx dy 2 ^ dx2 ôxôy dydx dy2 ôraknk p dfknkp

dTk drk

■ = 0;

nk j

pq, s _ p q, s ~

¡,j ¡,j , pq,s nk; ¡,j ~ nk; ¡,j ¡, j

+ nk; ¡,

p ¡, j - p ¡-i, j + pq,snjn±ii

hx

■ ¡, j

hx

pqj - pq:j-1 , pq,s nk;¡, j - nk;¡, j-1 _ D

+ nk; ¡,j " ' + p ¡,j

D

+-

2

nk ^ j

Pq, s pq,s

¡+1.j - p ¡-1, j + pq,s nk; ¡, j - nk; ¡-1, j

( pq, s pq, s A

n p ¡+1,j - p ¡,j , pq,s nk;¡+1,j~nk; ¡,j

h2 + p+1,j h2

■ ¡, 1

hx

(

nk-, ¡, j

Dq,s pq,s ^

p ¡,j ~ p ¡-1,j + pq,s nk; ¡,j ~nk;i-1,j

hxhy

• ¡, J

hxhy

D

+—

2

D

2

D

+ —

2

D 2

pq, s _pq, s

^¡,j-1 ^¡-1,j-1 Dq,s nk;¡, j-1 ~ nk;¡-1, j-1

nk; ¡,J-1 \h " + p1

( nq,s nq,s .. ^

hxhy

T>q, s Dq,s

p ¡, j ~ p ¡, j-1 , pq,s nk;¡, j - nk;¡, j-1

nk ■ ¡ j———+pqj

k j h h ¡, 1 ftyftx

y x

hyhx

nk v-1,1

pq, s pq, s

p + pq,s nk ; i-1, j - nk; ¡-1, j-1

hyhx

hyhx

pq,s _ pq,s

n ruj . pq,s n; ¡,J+1 ~ nk;¡,J

V nk; ¡, 1 + 1 hy2 ¡,1 + 1 hy2

D

nk ; ¡, j

pq,s _pq,s

p¡+1,j p¡,j—1 + pq,s nk;¡,j - nk;¡,j-1

Л

hy2

¡, i

hy2

+ nk ;¡, j p?,'

q, s j j

1 hT

pqs _ „q-1, s

q,s ¡,j ¡,j

+ nk;¡,j< j --

hT

s _ ^s-1

Л,1 Ai

- + nk;¡, jP^-¡, 1 "¡, 1 + n^j^

hT h

q,s P,j P,j = 0;

T

(12)

Ap PJ + Bp Pq,1s, j + CpPj + DpPqf, j + Ep pqjs+1 + Fp pqj1s, j _1 + Gp p«j-1,s + Hp J = 7p.

Численное моделирование горения аэровзвеси частиц алюминия... Здесь прогоночные коэффициенты имеют вид

Ap =

2nk,i,j - nk,i,j 2nk,i,j - nk,i-\,j 2nk,i,j - Hk,i,j-i

T hx hy

_ d 2nt,i,j - nk,i-i,j + nk,i-i,j _ d j - nk,i-i,j + nk,i+i,j

2hX hxhy

, 2nk,i,j - nkÂj_i + nki,j+i 2&f'j - aj 2fQj - fQ

D-^-+ nkij-+ nkij

2hy k,',j hT j hT

C Dnki,j Dnki,j D _ Dnk,i+i,j

Cp =" hy + 2hxhy " 2h2y , Dp = " 2h2 ;

E _ D jnk,i,j+i - nk,i,j ) F _ Dnk,i-i,j

Ep " 2h2 ' Fp " 2hxhy ;

q,s rq, s ^ q,s

G nk,i,j®j,j H nk,i,jfi j j nk,i,jPj,j

GP ~ 7 , HP ~ 7 , JP ~

hT hT x

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

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

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

О _ Tk - Tk0 . û _ rk - rk min . fi _ mi . fi _ x - x0 . 0 _ y - yo

~ rp T ' r ~ > Ъщ — . Ux — > Oy — ,

Tkmax _ Tk0 rkmax _ rkmin X» _ X0 ymax - y0

где rkmin = 0,5 мкм; rkmax = 6 мкм; Tkmax = 2300 K; X0 = 0; xœ = 5 см; y0 = 0. ymax = 3 см.

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

Результаты численного моделирования. Вид двумерной функции ПРВ во фронте пламени приведен на рис. 2, а. Зависимость характеризуется наличием двух локальных максимумов, соответствующих состоянию частиц, близкому к начальному, и горению в парофазном режиме, причем доля частиц составляет примерно 22 %. Наиболее вероятно нахождение частицы во фронте пламени с параметрами 9т е[0,9;1] и 9r е[0;0,07].

Вид двумерной функции ПРВ во фронте пламени в проекциях на оси 9т и 9y и оси 9r и 9x соответственно представлен на рис. 2, б, в. Частицы сначала распределены по нормальному закону, но потом под действием компонента скорости w перемещаются и прилипают к одной из стенок.

Рис. 2 (окончание). Вид функции ПРВ во фронте пламени в проекции на оси 6Г

и 6х (в)

Распределения суммарного объемного тепловыделения по пространственным координатам х и у при значениях коэффициента а = 0,5 и 1 приведены на рис. 3. При увеличении коэффициента избытка окислителя до 1 (стехиометри-ческий состав) снижается максимальное значение объемного тепловыделения О и уменьшается градиент всех параметров во фронте пламени.

q, 10" 8, Вт/см3

q, 10" s, Вт/смJ

Рис. 3. Распределения суммарного тепловыделения при коэффициентах избытка окислителя а = 0,5 (1) и а = 1 (2) в проекциях на оси 6х (а) и 6у (б)

Для верификации численной модели проведено сравнение полученных численных результатов математического моделирования с экспериментальными данными, представленными в работе [15]. Получено удовлетворительное качественное согласование по градиенту изменения температуры и тепловыделения в камере сгорания. Максимальное расхождение вычисленных и экспериментальных значений градиентов не превысило 9,7 %.

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

ЛИТЕРАТУРА

1. Исследование реакции окисления активированного алюминия водой — метод получения водорода / А.В. Пармузина, О.В. Кравченко, Е.И. Школьников и др. // Изв. РАН. Сер. Химия. 2009. № 3. С. 483-488.

2. Ягодников Д.А., ред. Актуальные проблемы ракетного двигателестроения. М.: Изд-во МГТУ им. Н.Э. Баумана, 2017. 295 с.

3. Милёхин Ю.М., Ключников А.Н., Бурский Г.В. Энергетика ракетных двигателей на твердом топливе. М.: Наука, 2013. 207 с.

4. Davis A. Solid propellants: the combustion of particles of metal ingredients // Combust. Flame. 1963. Vol. 7. No. 4. Р. 359-367. DOI: 10.1016/0010-2180(63)90212-8

5. Gnanaprakash K., Chakravarthy S.R., Sarathi R. Combustion mechanism of composite solid propellant sandwiches containing nano-aluminium // Combust. Flame. 2017. Vol. 182. P. 64-75. DOI: 10.1016/j.combustflame.2017.04.024

6. Yagodnikov D.A., Gusachenko E.I. Effect of an external electric field on the disperse composition of condensed products of aluminum particle combustion in air // Combust. Explos. Shock Waves. 2002. Vol. 38. Iss. 4. P. 449-455. DOI: 10.1023/A:1016215416889

7. Study of aluminum particle combustion in solid propellant plumes using digital in-line holography and imaging pyrometry / Chena Yi, D.R. Guildenbechera, K.N.G. Hoffmeistera, M.A. Coopera, et al. // Combust. Flame. 2017. Vol. 182. P. 225-237.

8. Щетинин Г.А. Численное моделирование горения аэровзвеси частиц алюминия с использованием функции плотности распределения вероятности // Политехнический молодежный журнал. 2017. № 11 (16). DOI: 10.18698/2541-8009-2017-11-204

9. Гавин А.Б., Медведев В.А., Наумов Н.А. Модель двухфазной турбулентной струи с учетом гетерогенного горения частиц // Физика горения и взрыва. 1988. Т. 24. № 3. С. 12-17.

10. Ягодников Д.А., Воронецкий А.В. Влияние скоростной неравновесности на особенности распространения ламинарного пламени в аэродисперсной среде // Физика горения и взрыва. 1992. Т. 28. № 5. С. 38-44.

11. Вилюнов В.Н., Ворожцов А.Б., Фещенко Ю.В. Моделирование двухфазного течения смеси газа с горящими частицами металла в полузамкнутом канале // Физика горения и взрыва. 1989. Т. 25. № 3. С. 39-43.

12. Басевич В.Я., Володин В.П., Перегудов Н.И. Определение функции плотности вероятности температуры расчетом турбулентного пламени по мгновенным параметрам // Физика горения и взрыва. 1990. Т. 26. № 6. С. 22-26.

13. Беляев Ю.В., Бриль А.П., Жданович О.Б., Ходыко Ю.В. О верификации модельных распределений вероятностей в турбулентных нагретых струях // Физика горения и взрыва. 1990. Т. 26. № 6. С. 92-97.

14. Ягодников Д.А. Статистическая модель горения боровоздушной смеси в турбулентном потоке // Физика горения и взрыва. 1996. Т. 32. № 6. С. 29-46.

15. Ягодников Д.А. Воспламенение и горение порошкообразных металлов. М.: Изд-во МГТУ им. Н.Э. Баумана, 2009. 432 с.

Романова Татьяна Николаевна — канд. физ.-мат. наук, доцент кафедры «Программное обеспечение ЭВМ и информационные технологии» МГТУ им. Н.Э. Баумана (Российская Федерация, 105005, Москва, 2-я Бауманская ул., д. 5, стр. 1).

Ягодников Дмитрий Алексеевич — д-р техн. наук, профессор, заведующий кафедрой «Ракетные двигатели» МГТУ им. Н.Э. Баумана (Российская Федерация, 105005, Москва, 2-я Бауманская ул., д. 5, стр. 1).

Щетинин Григорий Александрович — студент кафедры «Программное обеспечение ЭВМ и информационные технологии» МГТУ им. Н.Э. Баумана (Российская Федерация, 105005, Москва, 2-я Бауманская ул., д. 5, стр. 1).

Просьба ссылаться на эту статью следующим образом:

Романова Т.Н., Ягодников Д.А., Щетинин Г.А. Численное моделирование горения аэровзвеси частиц алюминия с использованием функции плотности распределения вероятности и двух пространственных координат // Вестник МГТУ им. Н.Э. Баумана. Сер. Естественные науки. 2018. № 6. C. 75-91. DOI: 10.18698/1812-3368-2018-6-75-91

NUMERICAL SIMULATION OF ALUMINUM PARTICLE AEROSOL COMBUSTION EMPLOYING A PROBABILITY DENSITY FUNCTION AND TWO SPATIAL COORDINATES

T.N. Romanova D.A. Yagodnikov G.A. Shchetinin

rtn@bmstu.ru daj@bmstu.ru gashetinin@yandex.ru

Bauman Moscow State Technical University, Moscow, Russian Federation

Abstract

The paper considers a non-steady-state two-dimensional model employing a probability density function describing temperature distribution and particle radii to investigate two-phase reactive flows formed from polydisperse aluminum particle aerosols. Our numerical computation technique uses six- and eight-diagonal matrices and the Gaussian method adapted to those. This modification aided in making the computations orders of magnitude faster. We obtained particle state distributions in the zones of maximum heat generation and aluminum burnout for the laminar aerosol flow mode. We detected specific features of the spatial flow structure characterised by development of reverse flow zones filled with burning aluminum particles

Keywords

Ignition, combustion, powdered aluminum, polydisperse aerosol, two-dimensional probability density function, mathematical model, numerical calculation

Received 09.04.2018 © BMSTU, 2018

REFERENCES

[1] Parmuzina A.V., Kravchenko O.V., Shkolnikov E.I., et al. Research on oxidizing reaction of activated aluminium by water — a method for obtaining hydrogen. Izvestiya RAN. Ser. Khimiya [Russian Chemical Bulletin], 2009, no. 3, pp. 483-488 (in Russ.).

[2] Yagodnikov D.A., ed. Aktual'nye problemy raketnogo dvigatelestroeniya [Actual problems of rocket engine technology]. Moscow, Bauman MSTU Publ., 2017. 295 p.

[3] Milekhin Yu.M., Klyuchnikov A.N., Burskiy G.V. Energetika raketnykh dvigateley na tverdom toplive [Energetics of solid-propellant rocket engine]. Moscow, Nauka Publ., 2013. 207 p.

[4] Davis A. Solid propellants: the combustion of particles of metal ingredients. Combust. Flame, 1963, vol. 7, no. 4, pp. 359-367. DOI: 10.1016/0010-2180(63)90212-8

[5] Gnanaprakash K., Chakravarthy S.R., Sarathi R. Combustion mechanism of composite solid propellant sandwiches containing nano-aluminium. Combust. Flame, 2017, vol. 182, pp. 64-75. DOI: 10.1016/j.combustflame.2017.04.024

[6] Yagodnikov D.A., Gusachenko E.I. Effect of an external electric field on the disperse composition of condensed products of aluminum particle combustion in air. Combust. Explos. Shock Waves, 2002, vol. 38, iss. 4, pp. 449-455. DOI: 10.1023/A:1016215416889

[7] Yi Chena, Guildenbechera D.R., Hoffmeistera K.N.G., Coopera M.A., Stauffachera H.L., Olivera M.S., Washburnb E.B. Study of aluminum particle combustion in solid propellant plumes using digital in-line holography and imaging pyrometry. Combust. Flame, 2017, vol. 182, pp. 225-237.

[8] Shchetinin G.A. Numerical simulation of aluminium particle aerosol combustion employing a probability density function. Politekhnicheskiy molodezhnyy zhurnal [Politechnical Student Journal], 2017, no. 11 (16) (in Russ.). DOI: 10.18698/2541-8009-2017-11-204

[9] Gavin A.B., Medvedev V.A., Naumov N.A. Model of a two-phase turbulent jet with heterogeneous particle combustion taken into account. Combust. Explos. Shock Waves, 1988, vol. 24, iss. 3, pp. 271-276. DOI: 10.1007/BF00750604

[10] Yagodnikov D.A., Voronetskii A.V. Effect of velocity nonequilibrium on the laminar flame propagation characteristics in an air-dispersed medium. Combust. Explos. Shock Waves, 1992, vol. 28, iss. 5, pp. 484-490. DOI: 10.1007/BF00755719

[11] Vilyunov V.N., Vorozhtsov A.B., Feshchenko Yu.V. Modeling of two-phase flow of a gas mixture with burning metal particles in a semienclosed channel. Combust. Explos. Shock Waves, 1989, vol. 25, iss. 3, pp. 296-300. DOI: 10.1007/BF00788801

[12] Basevich V.Y., Volodin V.P., Peregudov N.I. Determining the probability density function of the temperature by calculating a turbulent flame from the instantaneous parameters. Combust. Explos. Shock Waves, 1990, vol. 26, iss. 6, pp. 640-644. DOI: 10.1007/BF00786500

[13] Belyaev Yu.V., Bril' A.P., Zhdanovich O.B., Khodyko Yu.V. Verification of model probability distributions in turbulent heated jets. Combust. Explos. Shock Waves, 1990, vol. 26, iss. 6, pp. 709-714. DOI: 10.1007/BF00786512

[14] Yagodnikov D.A. Statistical model of flame-front propagation in a boron-air mixture. Combust. Explos. Shock Waves, 1996, vol. 32, iss. 6, pp. 623-636. DOI: 10.1007/BF02111563

[15] Yagodnikov D.A. Vosplamenenie i gorenie poroshkoobraznykh metallov [Inflammation and combustion of metal powder]. Moscow, Bauman MSTU Publ., 2009. 432 p.

Romanova T.N. — Cand. Sc. (Phys.-Math.), Assoc. Professor, Department of Computer Software and Information Technologies, Bauman Moscow State Technical University (2-ya Baumanskaya ul. 5, str. 1, Moscow, 105005 Russian Federation).

Yagodnikov D.A. — Dr. Sc. (Eng.), Professor, Head of Department of Rocket Engines, Bauman Moscow State Technical University (2-ya Baumanskaya ul. 5, str. 1, Moscow, 105005 Russian Federation).

Shchetinin G.A. — student, Department of Computer Software and Information Technology, Bauman Moscow State Technical University (2-ya Baumanskaya ul. 5, str. 1, Moscow, 105005 Russian Federation).

Please cite this article in English as:

Romanova T.N., Yagodnikov D.A., Shchetinin G.A. Numerical Simulation of Aluminum Particle Aerosol Combustion Employing a Probability Density Function and Two Spatial Coordinates. Vestn. Mosk. Gos. Tekh. Univ. im. N.E. Baumana, Estestv. Nauki [Herald of the Bauman Moscow State Tech. Univ., Nat. Sci.], 2018, no. 6, pp. 75-91 (in Russ.). DOI: 10.18698/1812-3368-2018-6-75-91