CC BY
37
УДК 517.948
ПОСТРОЕНИЕ И ЧИСЛЕННАЯ РЕАЛИЗАЦИЯ МАТЕМАТИЧЕСКОЙ МОДЕЛИ ПРОЦЕССА ОБРАЗОВАНИЯ ПАРОГАЗОВОЙ СМЕСИ
И.И. Валов, А,Д. Жаботпинский, Б.М. Тюлькин
A CONSTRUCTION AND A NUMERIC REALIZATION OF STEAM-GAS STREAM MATHEMATIC MODEL
I.I. Valov, A.D. Zhabotinsky, В.M. Tulkin
В работе построена математическая модель образования парогазовой смеси и предложен алгоритм численного метода, предназначенного для расчета одномерных парогазовых потоков. Предложен стохастический алгоритм, позволяющий учесть случайный характер формирования подаваемых в поток капель в соответствии с действующим законом распределения их по диаметрам.
Ключевые слова: образование парогазовой смеси, численный метод, случайное распределение.
Mathematic model of steam-gas stream formation has been constructed, and numeric method of one-dimention two-phases streams calculation has been suggested. Stochastic algorithm to account the drop diameter in accordance to the current random distribution law has been suggested.
Keywords: a steam-gas stream formation, a numeric method, a random distribution.
Продукты сгорания ракетных двигателей имеют достаточно высокую температуру, что приводит к необходимости тепловой защиты близрасположенных конструктивных элементов. В настоящей работе предлагается способ снижения температуры путем впрыска воды в дозвуковой поток горячих газов.
Смесь, образующаяся в результате впрыска жидкости в поток горячего газа, представляет собой двухфазную трехкомпонентную среду, состоящую из газа, капель жидкости и образующегося при ее испарении пара. Для каждой составляющей определены свои макроскопические параметры, присущие каждой фазе.
Построение инженерной математической модели образования и дозвукового движения парогазовой смеси проведем при следующих сделанных предположениях:
1) газ считается инертным, то есть не вступающим в химические реакции с парами впрыскиваемой жидкости. В процессе смесеобразования его состав не претерпевает существенных изменений;
2) перемешивание паров впрыскиваемой жидкости и инертного газа происходит достаточно быстро, так что выполняется условие теплового равновесия;
3) стенки газохода считаем теплоизолированными;
4) нагрев и испарение капель впрыскиваемой жидкости происходит только за счет тепла, отбираемого от газовой фазы;
5) занимаемый каплями объем мал. Вследствие этого пренебрегаем процессами столкновения капель;
6) процессами дробления капель также пренебрегаем;
7) впрыскиваемая жидкость предполагается несжимаемой;
8) скоростным запаздыванием капель пренебрегается;
9) размер капель жидкости много меньше расстояний, на которых макроскопические параметры меняются существенно;
10) размер капель много меньше расстояния между ними, поэтому члены, характеризующие межфазное взаимодействие, находятся исходя из данных о теплообмене и массообмене для одиночной капли;
11) инертный газ и пар считаются идеальными с постоянными удельными теплоемкостями.
На основе сделанных предположений система уравнений для описания движения парогазовой смеси имеет вид:
Здесь р - плотность; и - скорость; Р - давление; Е - полная внутренняя энергия единицы объема; е - внутренняя энергия единицы массы газа; Д - газовая постоянная. Индекс г относится к инертному газу; п к пару; гф к газовой фазе. Величины, не имеющие индекса, относятся ко всей смеси в целом. Под плотностью рв понимается масса капель воды, содержащихся в единице объема смеси.
Первые два уравнения в системе (1) представляют собой запись закона сохранения массы. Первое уравнение показывает, что в потоке нет источников инертного газа. Второе характеризует тот факт, что состав парогазовой смеси изменяется за счет испарения капель жидкости в потоке. Третье уравнение - запись закона сохранения импульса в предположении об отсутствии скоростного запаздывания капель. Четвертое показывает, что отбираемое от газовой фазы тепло расходуется на нагрев и испарение капель. Система (1) дополняется уравнением состояния (2).
Для того, чтобы система уравнений (1 - 2) была замкнутой, требуется дополнить ее выражениями, определяющими величины I ж <5-
Величины / и <3 характеризуют массообмен и теплообмен между фазами. Их определение основано на исследовании процесса испарения единичной капли в газовой среде.
Описанию процесса испарения единичной капли посвящено множество работ [1, 2, 3]. Воспользуемся результатами, полученными Д. А. Франк-Каменецким [3]. Для единичной капли приток массы и энергии определяется для двух основных случаев:
1. Температура капли меньше температуры кипения воды:
Теплота, отбираемая у газовой фазы, идет на нагрев капли. В предположении, что ввиду малости диаметра капли ее температура постоянна по всему объему, дифференциальное
(1)
Д(Ягф) I 4(£гф+-РЬ.
Р = (рКГ) гф = (ргКт + рпЛп)Тгф,
(2)
Ргф — Рг + Рп р — Ргф "Ь Рв>
(ре)гф = (ре)г + (ре) п, ег = (Су)г ■ Тг еп — (Су)п • Тп.
уравнение изменения температуры капли записывается в следующем виде:
ёТ 6А
« (3)
Секундный приток энергии к капле составляет:
д = 7г£>2А(Тгф-Т). (4)
Ввиду кратковременности процесса нагрева до температуры кипения пренебрежем возможным диффузионным притоком массы к газовой фазе:
3 = О- (5)
2. Температура капли равна температуре кипения воды: Теплота, отбираемая у газовой фазы, идет на испарение капли. Приток массы к газовой фазе составляет
• ~ г.А-ЛГп
3 = 27гД------/п
(Ср) п
(6)
при этом от газовой фазы отбирается количество тепла
Ч^-З-Ь. (7)
Также заметим, что следующие характеристики воды не являются константами:
- температура кипения воды Т5 при различных давлениях;
- удельная теплоемкость воды Св при различных температурах;
- удельная теплота парообразования Ь.
Значения этих характеристик определяются экспериментальным путем, результаты измерений представляются в виде таблиц зависимостей от параметра [4, 5]. На основании анализа представленных в этих таблицах данных были подобраны следующие аналитические зависимости:
В диапазоне изменения давлений от 1 до 100 атмосферных зависимость температуры кипения воды от давления Т$(Р) может быть аппроксимирована с помощью следующей формулы:
Г5(Р) = 274 + А ■ Рв,
(А = 5,796471166
[В = 0,2470900564,
что соответствует диапазону температур (373...584 К). Здесь Т$ выражается в Кельвинах, давление в Паскалях.
В диапазоне температур от 293 до 593 К (т. е. от 20 до 320 °С) удельная теплоемкость воды в зависимости от температуры СВ(Т) выражается как:
СВ(Т) = А(Т - 293)3 + В{Т - 293)2 + С(Т - 293) + А
'А = 0,0002 В = -0,057 (7 = 4,6015 Б = 4152,257,
удельная теплоемкость Св выражена в Дж/(кг*К), температура Т в Кельвинах.
В диапазоне температур от 373 до 573 К (от 100 до 300 °С) для удельной теплоты парообразования в зависимости от температуры насыщения Ь(Т) подобрано выражение
Удельная теплота парообразования Ь выражается в Дж/кг, температура насыщения Т в Кельвинах.
Система (1) - (2) с учетом выражений (3) - (7) представляет собой систему дифференциальных уравнений в частных производных, описывающую движение парогазокапельной смеси в одномерной постановке. Для ее решения предлагается применить численную схему Лакса-Вендроффа [6]. Для рассматриваемой задачи схема принимает следующий вид:
Шаг по времени Д£ является переменным и выбирается исходя из условия Куранта -Фридрихса - Леви:
Значения и и с рассматриваются во всех узлах сетки на предыдущем временном слое. Первый шаг - вычисление значений в промежуточных точках сетки:
Ь(Т) = А-Т2 + В Т + С,
'А = -10,90392 < В = 6059,66832 С = 15112363,502.
V 1
Аt < min
Здесь Cr - число Куранта, выбираемое в диапазоне (0... 1).
Аі
(Егф)Г+°о?5 = 5 [(-Егф)І + (ЗДІГ1]
Аі
2Ах
Аі
[<(№гф)п + К) - <ЧЕГфУп1 - Р,
[Я^ + ЯЇ]',
(л)ІГніл = 5 [(Л)ІГ+1 + (л)І+і] ;
□¿-НО,5 п+0,5
^г—0,5 п+0,5
-¿+0,5 _ / чг+0,5 . / \г+0,5 .
Рп+0,5 — 1Ргф;п+0,5 ^РгфУп+0,5>
г—0,5 __ / \г—и,о | / \ъ—и,о
Рп+0,5 — ІРгфіп+0,5 + 'Аф-'п+О^
г—0,5
\і—0,5
и,
и,
/ \г+0,5
г+0,5 _ \Р ‘ ^)п+0,5 . п+0,5 — ¿+0,5 ’
Рп+0,5
/ \г—0,5
г-0,5 _ (Р ' и)п+0,5 . п+0,5 г—0,5 '
Рп+0,5
/ \г+и,& __ / \г+0,5 / \г+и,о
ІРп/п+0,5 — \Ргф^+0,5 ~~ иМп+О.б
/ \г—0,5 __ / \г—0,5 / \г—0,5
1Рп;п+0,5 — \Ргф/п+0,5 \Рт)п+0,Ы
(Рг)Йо^г + (Рп)а°о^п
\г+0,5
ІЇ? ^*+0,5
І-С'гфіп+о.б
/ \г+0,5 (ип"+о,б)2
'.Ргф>'п+0,5 2
(Ргта(^)Г + (^)п+оГ5(С,у)п
/ г —0,5 \2
с 171 'нг—0,5 с „ л*—0,5 1.ип+о,5^
У.-д'гф;п+0,5 УРгф/п+0,5 2 МЙо?5(ЗДг + (Рп)Г+°Д(Си)п'
Второй шаг - вычисление значений в регулярных точках сетки:
(Рг)п+0,5-^г + (Рп)п+0,5-^п
г-0,5
АІ
г+0,5
(Рг)п+1 — (Рг)п ((Рг) ' гл)п+0,5 ((Рг) ' и)п+0,.
(Ргф)п+1 — (Ргф)
(Р'“)п+1 = (Р'м)
ДІ
п~А^
Аі Ах
((ргф) • «)п+0?5 - ((^ф) • «)п+0І5
г—0,5
+ Ді • 7^+15
/■л . 7/2чг+0,5 , рг+0,5 VР и )п+0,5 ^ п+0,5
/ 2\*—0,5 г>г—0,5
(Р ' и )п+0,5 °п+0,5
Ді
г+0,5
г+0,5 , г}2-4"0,5\
г—0,5
ч г—0,5 , пі—0,5 >
дали = (£гф)»-<+Ї5((&ф)Г+о,5+кш - <+й((Егф)г+х5+тар
Г1)] -
^'Яп+1)
(Рп)п+1 ~~ (Ргф)п+1 (Рг)п+1! Рп+1 = (Ргф)?г+1 + (Рв)п+1;
г (р-«)п+1.
ип+1 - ~г >
Рп+1
rpi _ (Дтф)п+1 ~ (Ргф)п+1^ **2^ *■«+!
(Рг)п+1 (Су)г + (Рп)иЛСу)и + = [(Pr)n+i + (Pn)n+i] - П+1;
((Рг)п+1 ■ + (Рп)п+1 ■ ^п)
п+1
гг =
'-'эт -L1
((Ргф)п+гУ
Третий шаг - сглаживание решения по традиционной схеме метода: ql, q2 - вспомогательные величины. Cm - коэффициент сглаживания (обычно принимается равным 0,4).
ql = «Vil + - К+1\ - сги+1; q2 = |<+1| + сгп+1 - |<-\| - сг~+\;
MÎ+i = MU i + [((ft)« - - (0*)»+i - (ft)t\)l92|] ;
(Ргф)п+1 ~ (Ргф)п+1 + Cm[((Ргф)п"+1 (Ргф)п+1)1?1| ((Ргф)п+1 (Ргф)п+1)к2|] ;
(р ■ и)гп+1 = (р • и)гп+1 + Cm^ [((р • и);+\ - (р ■ u)ln+l)\ql\ - ((р • и)гп+1 - (р • «)^+\)|д2|] ;
(ад+i = (ЗД+i + Ст^ [((Егфуп\\ - (^гф)п+1)|д1| - ((Ягфй+х - (ЕгфУп+\)№} ;
(Рп)п+1 = (Ргф)п+Х “ (Pr)n+lî Рп+1 (ргф)тг+1 "I- (рв) П+1 >
г _ (Р • ц)п+1 ,
un+l i )
Рп+1
i 12
P?+l = [{Pvt+xR + {Pu)UiR] -П+ъ
г _ /((Рг)п-ц ' кг + (Рп)„+х ' &п) ' Рп+1
Сп+1 “ V ((Ргф)^г)2
Как непосредственно следует из рассмотрения схемы численного метода, для ее реализации на каждом следующем расчетном шаге по времени требуется определить величины (рв)п+1, ■1п+1 и <&+!•Их нахождение основывается на исследовании процесса испарения единичной капли в газовой среде и для ячейки с индексом (г,п) эти величины представляют собой:
(Рв)£г+1 ~ сумма масс всех капель, заключенных в ячейке расчетной сетки, деленная на объем этой ячейки;
Гп+1 - сумма притоков массы от всех капель, заключенных в ячейке расчетной сетки, деленная на объем этой ячейки;
Яп+1 сумма притоков энергии от всех капель, заключенных в ячейке расчетной сетки, деленная на объем этой ячейки.
При программной реализации метода для нахождения {рв)гп+1> ^п+г и Оп+1 предлагается ввести в рассмотрение новый объект - «Капля », обладающий следующими основными характеристиками:
- положение капли (т. е. ее координата ж);
- данные о текущих диаметре (и как следствие массе) и температуре капли;
- данные о текущем притоке массы и текущем притоке энергии от капли, определяемые в соответствии с (4) - (7).
Все участвующие в расчете капли помещаются в динамический список. Первоначальная температура каждой капли принимается равной температуре воды в емкости.
Формирование капель происходит следующим образом:
Вводится в рассмотрение некоторая величина - «ожидаемая масса» М. На каждом шаге по времени Д£ определяется текущий расход каждой форсунки гп и «текущая мас-са»(Мт)г+1 = {Mт)i + ш •
Согласно теории форсунок [7], расход каждой форсунки т определяется как:
т
М ■ >5ф\/2 • А-Рф • рж,
¡j, < 1 - коэффициент расхода форсунки, 5ф - площадь проходного сечения форсунки, ДРф
- перепад давлений на форсунке, рж - плотность впрыскиваемой через форсунку жидкости.
Когда текущая масса превышает ожидаемую, капля считается сформированной и вводится в расчет (добавляется в динамический список). «Текущая масса» обнуляется и назначается новая «ожидаемая масса».
Поскольку диаметр капель жидкости, подаваемой в поток горячего газа через форсунку, не является постоянным, а носит случайный характер, при программной реализации метода требуется дополнительно реализовать алгоритм, возвращающий случайное значение ожидаемого диаметра капли в соответствии с действующим законом распределения.
Как показывают проведенные опыты [7], при впрыске жидкости в газовый поток распределение впрыскиваемых капель по их диаметрам в факеле распыла форсунки подчинено нормальному закону распределения. Плотность нормального распределения выражается формулой
1 AD-Md?
ф{0) = 775е -2^—
При этом математическое ожидание М,/ и среднее квадратичное отклонение а диаметра капли зависят от перепада давлений на форсунке и ее конструктивных особенностей, а также характеристик газа в потоке и подаваемой жидкости.
Одна из теорий распада струи жидкости в газовом потоке гласит, что распад струи происходит по причине воздействия внешних сил на ее поверхность. При движении струи происходит трение ее поверхности о сжатые газы. В результате действия сил трения и сил поверхностного натяжения струя распадается на отдельные капли. При этом математическое ожидание диаметра капли Мд определяется по полуэмпирической формуле [7]:
Md = 3.33
« Рж
Д-Рф Ргф (р Ргф)
Здесь а - коэффициент поверхностного натяжения подаваемой жидкости.
Среднее квадратичное отклонение диаметра капли а определим как а = кв • кй -коэффициент вариации. Примем этот коэффициент постоянным (то есть не изменяющимся в процессе всего расчета) и зависящим только от конструктивных особенностей форсунки. Пусть кв находится в пределах (0... |) (в силу того, что диаметр капли П не может принимать отрицательные значения).
Чтобы отразить в математической модели случайный характер диаметра образующейся капли, следует реализовать алгоритм, возвращающий ожидаемое значение диаметра, исходя из известных значений его математического ожидания и среднеквадратичного отклонения.
Исходя из свойств нормального распределения и воспользовавшись «правилом 3 сигма», будем предполагать, что диаметр капли может варьироваться в пределах [М({ — За... М(/ + За] (и не выходит из этого диапазона).
Зададимся неким достаточно большим целым числом N. Разобьем отрезок [Ма — 3 а... М(1 + За] на п <С N равных элементарных отрезков, к — ^ - длина элементарного отрезка. Каждому элементарному отрезку с индексом г поставим в соответствие три целых числа: К{, (К тт)г, (К тах)^ К% выберем таким образом, чтобы
Ki
Md—3o-+ih
/ *{x)dx■
N
Md—3cr+(i—l)h
В соответствии с формулой трапеций имеем следующее выражение:
ф{Ма - За + (г - 1 )h) + ф(Ма - За + ih)
Ki =
h-N
В данном случае ||-|| обозначает взятие целой части числа.
Величины (К min)j, (К max),; определим как:
(К min)o = О, (К тах)о = Ко,
(К min)j = (К max)j_i i = 1,2,... n,
(К max)j = (К min)j + Ki i = 1,2,... n.
(Это же выражение задает и последовательность их вычисления).
На следующем шаге алгоритма сгенерируем случайное число R, изменяющееся в пределах [0.,.(iCmax)n]. При реализации программы на ЭВМ можно воспользоваться стандартной для многих языков программирования функцией random или одним из методов генерации псевдослучайных чисел.
Найдем индекс элементарного отрезка, для которого (К пііп)г < Я < (К тах)г. Возвращаемый диаметр В определим как координату середины этого отрезка:
В = Ма - За + (і - 0 • к.
В Таблице представлено сравнение расчетных данных с экспериментальными (проводится сравнение температурных режимов в рассматриваемом объеме).
Таблица
Температура в изучаемом объеме
Время, с Температура расчетная, К Температура измеренная, К
0 297 297
0,05 587 510
од 605 548
0,15 617 556
0,2 624 560
0,25 620 563
0,3 617 561
0,35 610 562
0,4 622 559
0,45 625 564
Литература
1. Нигматулин, Р.И. Основы механики гетерогенных сред / Р.И. Нигматулин. - М. : Наука, 1978.
2. Фукс, H.A. Испарение и рост капель в газообразной среде / H.A. Фукс. - М. : Изд-во Акад. наук СССР, 1958.
3. Франк-Каменецкий, Д.А. Диффузия и теплопередача в химической кинетике / Д.А. Франк-Каменецкий. - М. : Наука, 1987.
4. Кухлинг, X. Справочник по физике / X. Кухлинг. - М. : Мир, 1982.
5. Вукалович, М.П. Таблицы теплофизических свойств воды и водяного пара / М.П. Ву-калович. - М. : Изд-во стандартов, 1969.
6. Тюлькин, Б.М. Расчет ударно-волновых и струйных течений газа численным методом типа Лакса - Вендроффа / Б.М. Тюлькин // Сборник материалов XXIX студенческой научной конференции «Студент и научно-технический прогресс». - 2005. - С. 18 - 33.
7. Алемасов, В.Е. Теория ракетных двигателей / В.Е. Алемасов. - М.: Машиностроение, 1989.
Валов Илья Игоревич, кафедра прикладной математики, Математический факультет, Челябинский государственный университет,vii_1982@mail.ru.
Жаботинский Анатолий Данилович, профессор, доктор технических наук, кафедра летательных аппаратов и автоматических установок, Аэрокосмический факультет, Южно-Уральский государственный университет, src@makeyev.ru.
Тюлькин Борис Михайлович, профессор, доктор технических наук, кафедра прикладной математики, Математический факультет, Челябинский государственный университет.
Поступила в редакцию 21 июля 2010 г.
