Газовая и волновая динамика - часть 32

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

Контактные разрывы, как зоны больших градиентов плотности, хорошо видны

на шлирен-фотографиях рис. 3 – 4. Горячие точки перед фронтом пламени, порож-
дающие волны детонации (рис. 3) или дефлаграции (рис. 4a), образуются именно
на контактных поверхностях, являющихся результатом взаимодействия догоняющих
друг друга головных ударных волн. Теоретическое объяснение этого факта было
впервые предложено в работе (Смирнов, Панфилов, 1995), в которой было показано,
что само воспламенение в области ударно сжатого газа происходит в отдельном эк-
зотермическом центре, характеризуемом минимальным значением времени задержки
воспламенения. Последнее зависит от температуры смеси и от физического време-
ни экспозиции смеси при данной температуре. Температура газа перед контактной
поверхностью, сжатого головной ударной волной, выше, чем температура газа за кон-
тактной поверхностью, прошедшего последовательно две стадии ударного сжатия и
адиабатического расширения. Время индукции воспламенения с повышением темпе-
ратуры сокращается. В области течения между контактной поверхностью и ударной
волной слои газа, примыкающие к контактной поверхности, характеризуются самым
коротким временем индукции, так как они дольше всего находились при повышенной
температуре. Поэтому самовоспламенение происходит в горячей точке на контактной
поверхности. Эти рассуждения подтверждаются и результатами прямого численного
моделирования (Смирнов, Панфилов, 1995; Смирнов и др., 1999).

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

пературы и концентрации реагентов, установившиеся в окрестности зоны воспла-
менения, определяют, разовьется ли из данного очага воспламенения волна детона-
ции или дефлаграции (Зельдович и др., 1970). При возникновении волны детонации
реализуется картина течения, изображенная на рис. 3. При возникновении волны
дефлаграции, которая распространяется от места воспламенения с гораздо меньшей
скоростью (рис. 4a), возможно последующее воспламенение соседних «горячих то-
чек», до которых волна дефлаграции еще не успела дойти (рис. 4b).

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

Численные исследования процессов ПГД проводились с помощью системы уравне-
ний, полученной осреднением по Фавру системы уравнений движения многоком-
понентной смеси. Использовалась модифицированная ка-эпсилон модель турбулент-
ности. Для моделирования флуктуаций температуры к основным уравнениям ка-
эпсилон модели было добавлено третье уравнение, описывающее динамику средне-
квадратичного отклонения температуры (Смирнов и др., 2001). Члены, отвечающие
за производство и диссипацию этого параметра, а также нелинейные по температуре
члены, отвечающие за скорость химических реакций, моделировались с использова-

132

нием метода квадратур Гаусса.

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

параметров имеет следующий вид:

∂

t

(ρ) +

∇ · (ρu) = 0,

(1)

∂

t

(ρY

k

) +

∇ · (ρuY

k

) =

−∇ · I

k

+ ˙ω

k

,

(2)

∂

t

(ρ

u) +

∇ · (ρu ⊗ u) = ρg − ∇p + ∇ · τ,

(3)

∂

t

(ρE) +

∇ · (ρuE) = ρu · g − ∇ · pu − ∇ · I

q

+

∇ · (τ · u).

(4)

Уравнения (1) – (4) включают уравнение неразрывности для смеси в целом,

уравнения баланса массы для k-го компонента, уравнения количества движения
и энергии для смеси. Наряду с этими уравнениями должны также выполнять-
ся условия согласования:

k

Y

k

= 1,

k



I

k

= 0,

k

˙ω

k

= 0. Термическое и ка-

лорическое уравнения состояния для смеси газов имеют соответственно следу-
ющий вид: p = R

g

ρT

k

Y

k

/W

k

, E =

k

Y

k

(c

vk

T + h

0k

) + 

u

2

/2 + k. Турбулент-

ный тепловой поток в уравнении (4) представляется в виде суммы двух членов:



I

q

= 

J

q

+

k



I

k

(c

pk

T + h

0k

), где 

J

q

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

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

t

в рамках

ка-эпсилон модели определяется по следующей формуле: ν

t

= C

µ

k

2

/ε. Для опреде-

ления турбулентных потоков используются следующие соотношения:

τ = (µ + ρν

t

)(

∇u + ∇u

T

− (2/3)(∇ · u)U) − (2/3)ρkU,

(5)



I

k

=

−ρ(D + (ν

t

/σ

d

))

∇ · Y

k

,



J

q

=

−(λ +

k

c

pk

ρ

k

(ν

t

/σ

t

))

∇ · T,

(6)

Рождение k-го компонента в результате химической реакции происходит со ско-

ростью ˙ω

k

, которая является суммой массовых скоростей ω

kj

всех реакций, протека-

ющих в газовой фазе. Члены ˙ω

k

в уравнениях, отвечающие за химические превраще-

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



T



. Для замыка-

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

133

энергии турбулентных пульсаций k, и диссипации ε, к которым добавляется третье
уравнение для определения среднеквадратичного отклонения температуры θ = T



T



.

∂

t

(ρk) +

∇ · (ρuk) = ∇ · ((µ + ρ(ν

t

/σ

k

))

∇k) + τ

t

:

∇u − ρε,

(7)

∂

t

(ρε) +

∇ · (ρuε) = ∇ · ((µ + ρ(ν

t

/σ

ε

))

∇ε) + (ε/k)(C

1ε

τ

t

:

∇u − C

2ε

ρε),

(8)

∂

t

(ρˆ

c

p

θ) +

∇ · (ρuˆc

p

θ) =

∇ · ((λ +

k

Y

k

c

pk

ρ(ν

t

/σ

k

))

∇θ) + P

θ

+ W

θ

− D

θ

,

(9)

где ˆ

c

p

=

k

c

pk

Y

k

члены, характеризующие производство P

θ

, W

θ

и диссипацию D

θ

величины θ, определяются следующими формулами:

P

θ

= 2ρ

k

Y

k

c

pk



ν

t

/σ

k



(

∇T )

2

,

W

θ

=

−

k

˙ω



k

T



h

0k

(10)

D

θ

= C

g

ρ

k

Y

k

c

pk

ε

k

θ

θ

m

− θ

.

При определении членов производства при химических реакциях W

θ

использо-

вался закон Аррениуса. Для расчета конкретных членов ˙ω



k

T



применялась техника

квадратур Гаусса:

T



A(T ) = θ

A(T +

√

3θ)

− A(T −

√

3θ)

2

√

3θ

.

(11)

Вид функции D

θ

, характеризующей диссипацию, определялся исходя из предпо-

ложения, что флуктуации температуры ограничены таким образом, чтобы величина
T = ¯

T + T



не могла принимать отрицательные значения. Поэтому в формулу введен

множитель 1/(θ

m

− θ). Чтобы оценить величину θ

m

примем во внимание следующие

соображения. Во-первых, для нормальной функции распределения вероятность того,
что среднеквадратичное отклонение вдвое превысит среднее значение, оказывает-
ся менее 1%. Во-вторых, среднее отклонение температуры в экспериментах (Филип,
1991) не превосходило половины максимума осредненного значения. Поэтому, для
θ

m

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

θ

m

= T

2



4.

Константы ка-эпсилон модели определялись стандартным образом

134

C

µ

= 0.09, C

1ε

= 1.45, C

2ε

= 1.92,

σ

d

= 1, σ

t

= 0.9, σ

k

= 1, σ

ε

= 1.3.

Константа C

g

в уравнении (10) для вычисления диссипации определялась с ис-

пользованием экспериментальных данных (Филип, 1991) следующим образом:

C

g

= 2.8.

Предполагалось, что газовая смесь содержит следующие компоненты: O

2

, C

n

H

m

,

CO, CO

2

, H

2

, H

2

O, N

2

. Химический потенциал топлива и его состав (n и m) могут

рассматриваться, как параметры модели. Учитывались следующие модельные брутто-
реакции:

C

n

H

m

+



n

2

+ σ

m

4

O

2

→ nCO + σ

m

2

H

2

O + (1

− σ)

m

2

H

2

,

CO +

1
2

O

2

→ CO

2

,

CO

2

+ M

→ CO +

1
2

O

2

+ M,

H

2

+

1
2

O

2

→ H

2

O,

H

2

O + M

→ H

2

+

1
2

O

2

+ M.

Здесь σ – доля воды при разложении углеводорода (зависит от типа топлива).

Число реакций будем обозначать K. Скорость образования k-го компонента при
одновременном протекании всех реакций определяется:

˙ω

k

=

B

j

=1

˙ω

kj

,

(12)

где ˙ω

kj

– скорость образования k-го компонента при j-ой реакции, B – число хими-

ческих реакций.

Тогда источниковый член в уравнении (9) для среднеквадратичной пульсации

температуры, определяемый химическими реакциями (10), примет вид:

W

θ

=

−T



B

j

=1

K

k

=1

h

0k

˙ω

kj

.

(13)

Каждый член ˙ω

kj

содержит функцию Аррениуса для определения скорости j-ой

реакции. Предполагается, что эта функция имеет следующий вид:

A

j

(T ) =



K

j

exp



−

T

aj

T

при T

≥ T

mj

,

0 при T < T

mj

,

(14)

где K

j

– предэкспоненциальный множитель, T

aj

– температура активации, T

mj

– ми-

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

135