CC BY
10
УДК 681.326
К вопросу об усилении и численной реализации полевой модели развития пожара
Ясинский Ф.Н., д-р физ.-мат. наук, Потемкина О.В., канд. хим. наук,
Сидоров С. Г., Ясинский И.Ф., кандидаты техн. наук,
Предложена методика вычислений для полевой модели пожара. Сформулированы способы, позволяющие уменьшить затраты компьютерного времени и увеличить точность результатов.
Ключевые слова: вязкость газовой смеси, дифференциальные уравнения, двухпараметрическая К-Е модель, моделирование горения.
On Strengthening and Numerical Realization of Developing Field Model of Fire
F.N. Yasinskiy, Doctor of Physics and Mathematics, O.V. Potemkina, Candidate of Chemistry,
S.G. Sidorov, I.F. Yasinskiy, Candidates of Engineering
The calculation method is offered for the fire spreading field model. The methods, allowing to decrease the compute time and increase the results exactness are formulated.
Key words: parallel programming, simulation, turbulence, sweep method.
Полевая модель пожара, согласно [1], опирается на следующую систему уравнений:
dp
yAu. ) = о,
dt ті дх Л '
д_
dt
j=1 j
з
и )+y dx- Hu )=-dp *dXj '<+pg',
Tii =v
(duj dujЛ
dXj дх,
з
- 2 цу dUL ,8j ,
3 k= dXk *
ЦрЬ ) + ydLipuJh ) = dL + y^_
dt j=1 dXj dt T^dx-
dp
d j=iJxj
d dh Cp dxj
(1)
(2)
(3)
R
dqR
dx
N
h = ho + j CpdT + yYkHk,
k=1
k=1
pk> 3
1
(4)
(5)
(6)
I (pYk ) + У ^ (Yk )-ydr. [pD %\ + SkAV
j=1 dxj
p -pRoT y Mr •
k=1 Mk
(8)
Здесь t- время; и, / = 1, 2, 3 - три компоненты скорости; р, Р, Л - плотность, давление и теплосодержание газовой смеси; Ук - весовые доли ее компонентов с молекулярным весом Мк и теплотой образования Нк, к = 1, 2, ..., М;
,у = 1, 2, 3 - лучистый поток в направлении
у; ду - ускорение свободного падения; Т - температура.
Все указанные выше дифференциальные уравнения могут быть приведены к единой форме:
dQ dQ d
-----+ У u,----------= y ---------
dt — j dx, — dx ,■
v dQ
Oq dxj
\
Q-
(9)
у = 1 ^ у у = 1^у
Здесь О любая из следующих величин:
О = (и1, и2, из, К, Т, Е, Р, р1,., рМ
где р1,., рМ - плотности компонент газовой
смеси.
Исключение представляет лишь уравнение неразрывности (1). Предлагаем заменить его на уравнение
( \
dp
dt
3 а 3
■'Уdt(puj)=У У
j=1 dxj
_d_
dxj
dp
op dx■.
[ p J \
(10)
1=1 ^1
Здесь стр - некоторая специально подобранная константа, учитывающая перенос массы в потоке за счет турбулентных пульсаций, а не только конвективно, как это делается обычно. Впрочем, замена (1) на (10) уже давно делается в теории химических лазеров [2].
Переход от (1) к (10) не только оправдан с физической точки зрения, но (что особенно важно) ведет к большим вычислительным упрощениям. Эти упрощения состоят не только в вытекающем в результате этого единообразии в определение всех искомых величин по уравнениям типа (10), но и в более быстрой сходимости.
В [1] ничего не говорится о предпочтительном методе вычисления турбулентной вязкости. Не говорится также и о расчете химической кинетики, и зависимости теплоемкостей компонентов от температуры. Опираясь на наш вычислительный опыт, рекомендуем поступить следующим образом.
Полная эффективная вязкость газовой смеси состоит из молекулярной vm0| и турбулентной vturb :
^ = р(чт01 + vtuгd ).
Для вычисления vturb наиболее удобны два метода: А.Н. Секундова и К-Е метод.
о
В модели А.Н. Секундова поле турбулентной вязкости описывается дифференциальным уравнением
3
д_
ді
'/іиіЬ
I и
і=і
д
дХі
ЧигЬ
з д | д
= 1 (''то! + Хх’игЬ))УигЬ
1=1дХ і К дхі
+ушЬ ((ь/8уто/ )й - У^піп (тоІ + вуїигЬ )УШЛ>. Здесь Г(•) - следующая функция:
Г (г ) = 0,2 г 22+1'47 •г + °'2;
у ’ г2 -1,47 • г +1
(11)
(12)
X, р, у - эмпирические константы (X = 2; Р = 0,06; у = 50).
На неподвижных твердых поверхностях Г ставится граничное условие |г = 0 . Чтобы
поле турбулентной вязкости сформировалось, нужно задать некоторое его начальное значение.
Еще более совершенной моделью для вычисления поля турбулентной вязкости является двухпараметрическая К-Е модель. Здесь К - кинетическая энергия турбулентных пульсаций, а Е - скорость диссипации кинетической энергии этих пульсаций. Модель задается следующими выражениями:
"іигЬ
= С,,
к_
Е
дК ^ дК ^ д
— + / и,-----------------= І -----------
ді дх • дх •
01 ,=1 ох] ,=1 ох]
увГГ
дК
дх,
иЕ V-1 дЕ V-1 д
— + І и,--------------= І ---------
ді дх • дх •
01 ]=1 дх] ]=1 0х]
Бк = V • й - Е;
Бе =(й - с2Е)• Е / К;
д
'^вгг дЕ
а к дх ]
(13)
(14)
(15)
(16) (17)
где Veгг =V, + Vturb; С,= 0,°9; С1 =1,44; С2 =192;
ак = 1; стє = -І,3.
Начальные и граничные условия обычно берут в виде
Ко = Ск І и?; Ео = се- Ко3/2, (18)
і
где ск, се - некоторые малые положительные постоянные.
Граничные условия на твердых поверхностях имеют вид
дК_
ду
= о; дЕ
г ду
= 0.
(19)
Здесь у - нормаль к неподвижной твердой поверхности.
Очевидно, что обе модели предлагают уравнения, удовлетворяющие универсальной форме (10).
Указанная выше система уравнений, дополняется уравнениями химической кинетики. Если Ха - символ компонента номер а, а ааг и Ьаг -стехнометрические коэффициенты г-реакции
І
аагХа
I Ьа
Xа
(20)
то изменение плотности ра в заданном жидком микрообъеме компонента номер а за счет реакций описывается системой
=
^ра
а
Ма ( Ьаг ааг )
d <вг
\ааг
- к.
Ьг
П
а У
d <вг
М,
\Ьг
(21)
(22)
а У
\Гг, кЬг - коэффициенты скоростей для
где К,г, Кь
прямой и обратной реакции, очень сильно зависящие от температуры:
-пг 1 а ,т\ (23)
,0
Кг = Кг0Гепхгр (-Ег / т),
где Кг, Ег, пг - некоторые постоянные, характерные для г-реакции.
Выделение или поглощение тепла на множестве реакций и соответствующее изменение температуры определяется источниковым членом для температурного уравнения 1 d юг
рт = сЦчг^т •
3 пт 3 „ ( ^
(24)
X д Т дх.
+Рт .
(25)
дТ дТ ^ д
— + > и,-------= > --
д^ 1 дх. ¿—‘ду-
дt 1=1 дх1 1=1 дх1 ^р-'7
Наша вычислительная практика позволяет рекомендовать следующую методику интегрирования описанной системы уравнений.
1. Расщепление универсального уравнения (10) для всех искомых величин на множество одномерных уравнений:
дО дО
— + иі —
ді 1 дхі
і = 1, 2, 3,
_д_
дхі
дО
дхі
dt
О
(26)
(27)
Локально-одномерное расщепление хорошо описано в [3].
В пределах малого шага по времени т эти уравнения интегрируются последовательно. Результаты, полученные при решении предыдущего уравнения, являются исходными для следующего. Интегрируются уравнения переноса методом скалярных прогонок по прямым, параллельным соответствующей оси. При построении разностных схем используются «противоточные производные», что существенно увеличивает устойчивость вычислений. Разностная схема для направления, например, х1 будет иметь вид
г
Цок « - ОК) +
+ 2Ъ ((( -))(°° - °К+1) +
+ ( + |<4|)(С - ОК-Ц,))- (28)
- ^ ((«¡+1,+*1,(1, - ОК *' )--(+V/-1,)*’ - ой;))■
Нижние индексы соответствуют номерам узлов трехмерной сетки, а верхний - номеру момента времени.
Это выражение приводится к виду, при
_ к +1
котором искомые щ, соответствующие новому моменту времени, собраны в левой части равенства
дКоК+1 + вКоК+1 + ГКпК+1 = пК (29)
Л//0/ -1/ + в/// О/1/ + С/11О/+11 = и/]1, (29)
где ЛК,вК,Сщ,Оу - коэффициенты, получающиеся при приведении (28) к форме (29).
Согласно правилам прогонки, решение ищется в виде
О-+ = Щ,С$+1 + МК. (30)
Прямая прогонка выполняется по формулам
К,, =-СК / ( + Л^), (31)
М,К+1У/ = ( -лЦм*)/( + ). (32)
А обратная - согласно (30).
После окончания операций с уравнениями переноса делается шаг по времени для уравнений химической кинетики (21), (22). Обычно они отличаются высокой жесткостью, мерой чего является отношение абсолютных значений собственных чисел (максимального к минимальному) якобиана правых частей кинетических уравнений.
Наилучшим в этом случае оказывается метод Гира [3].
При решении уравнений химической кинетики удобно воспользоваться асинхронным интегрированием, состоящим в том, что для
различных кинетических уравнений выбираются различные шаги по времени. Эти шаги являются кратными друг другу, что позволяет затем их вместить в некоторый объединяющий шаг. За счет такого подхода исключается ситуация, когда небольшая жесткая подсистема навязывает малый шаг по времени всей системе. При применении асинхронной технологии для значений переменных, соответствующих некоторому локальному времени, применяется интерполяция по Лагранжу. Отметим, что такая технология успешно применяется в молекулярной динамике [4] и в расчетах лазерных систем [5].
Главная трудность при моделировании горения состоит в большом числе уравнений химической кинетики и уравнений переноса для каждой из многочисленных компонент. Весьма обнадеживающим является следующее: возникшее в последнее время направление по созданию суперкомпьютеров на базе графических процессоров и программной системы CUDA обладает тем большей эффективностью, чем на большее число потоков удается разбить вычислительный процесс. Число таких потоков может быть очень велико. Это вселяет определенные надежды на прогресс в вопросе математического моделирования горения.
Список литературы
1. МЧС России, Приказ №382 от 30.06.2009. «Об утверждении методики определения расчетных величин пожарного риска в зданиях, сооружениях и строениях различных классов функциональной пожарной опасности».
2. Риврд У., Батлер Т., Фармер О. Нестационарные турбулентные течения химически реагирующих газовых смесей / Численное решение задач гидромеханики. - М.: Мир, 1977. - С. 184-193.
3. Калиткин Н.Н. Численные методы. - М.: Наука, 1978. - С. 512.
4. Street W.B., Tildesley, Multiple time-step methods in molecular dynamics // Molecular Physics. - 1978. -Vol. 35. - №03. - Р. 639-648.
5. О построении быстрых алгоритмов для задач релаксационной газовой динамики / Э.Л. Комаров, И.Н. Макаров, Н.Н. Романова, Ф.Н. Ясинский // Сб. «Физико-химическая кинетика в газовой динамике» / Институт механики МГУ. - М.: Изд-во МГУ, 1986. - С. 109-114.
Работа выполнена при поддержке Министерства образования и науки РФ ГК № 13 G25.31.0077.
Ясинский Федор Николаевич,
ГОУВПО «Ивановский государственный энергетический университет имени В.И. Ленина»,
доктор физико-математических наук, профессор кафедры высокопроизводительных вычислительных систем, телефон (4932) 26-98-29.
Потемкина Ольга Владимировна,
Ивановский институт государственной противопожарной службы МЧС России, заместитель начальника, телефон (4932) 26-16-24.
Сидоров Сергей Георгиевич,
ГОУВПО «Ивановский государственный энергетический университет имени В.И. Ленина»,
кандидат технических наук, доцент, зав. кафедрой высокопроизводительных вычислительных систем, телефон (4932) 26-98-29.
Ясинский Игорь Федорович,
Ивановская государственная текстильная академия,
кандидат технических наук, доцент кафедры прикладной математики и информационных технологий, e-mail: [email protected]
