Анализ термически возбуждаемых пульсационных колебаний газа в цилиндрической трубе Текст научной статьи по специальности «Физика»

АНАЛИЗ ТЕРМИЧЕСКИ ВОЗБУЖДАЕМЫХ ПУЛЬСАЦИОННЫХ КОЛЕБАНИЙ ГАЗА В ЦИЛИНДРИЧЕСКОЙ ТРУБЕ

Ж.М. САХАБУТДИНОВ *, О.С. КОЧНЕВА *, Г.И. ПАВЛОВ **

* Казанский государственный энергетический университет, г. Казань,

** филиал Военного артиллерийского университета, г. Казань

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

1. Введение

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

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

Условно среди этих задач можно выделить два подхода к изучению явления. В первом подходе параметры в «холодной» и «горячей» зонах удовлетворяют модели одномерного, идеального, нетеплопроводного газа. В зоне теплоподвода параметры газа меняются скачкообразно, решения в отдельных зонах «сшиваются» на основе специально выведенных уравнений. Как правило, при таком подходе уравнение энергии интегрируется отдельно от других уравнений. Такой подход наиболее полно отражен в работах [1, 2].

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

2. Основные дифференциальные уравнения сохранения

© Ж.М. Сахабутдинов, О.С. Кочнева, Г.И. Павлов Проблемы энергетики, 2004, № 3-4

Используется модель одномерного течения. Концы трубки остаются открытыми, давление на входе и выходе полагаются постоянными. В некотором сечении трубы имеется теплоподвод. Уравнения сохранения массы, импульса, тепловой энергии и уравнение состояния для идеального нетеплопроводного газа имеют вид [4]

др др ди

-------+ и — = —р — ,

д£ дх дх

(1)

ди дил

р| — + и— п Ы дх

.дР

дх

(2)

де дел

р| — + и —

1 дt дх у

ди

= — р ^ + Q' дх

(3)

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

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

р = рЯТ , (а) е = сиТ, (Ь) (4)

где я = си(у —1) - газовая постоянная; си - удельная теплоемкость при постоянном объеме; у - показатель адиабаты. Если ввести в рассмотрение удельную энтропию идеального газа

я = си 1п р — с р 1п р, то уравнение тепловой энергии (3) можно переписать в виде

рТ

(дэ дя Л

-----+ и —

Ы дх

= Q

Из этого представления видно, что энтропия элементарной частицы меняется только в случае неравенства нулю Q. Там, где отсутствует подвод теплоты, энтропия сохраняется, температуры соседних объемов не выравниваются из-за отсутствия вязкости и теплопроводности газа. Особенностью системы (1)-(4) в случае Q(x,£) = 0 является то, что первые два уравнения могут быть проинтегрированы отдельно от третьего. Физический смысл отделения третьего уравнения от первых двух сводится к тому, что акустические волны давления и скорости распространяются независимо от

«тепловых» волн энтропии. Это обстоятельство, как отмечалось во введении, накладывает свой отпечаток на метод решения задачи об устойчивости горения.

Используя второе уравнение (4), запишем уравнение тепловой энергии через температуру

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

Переходя в системе уравнений (1), (2), (3а) и (4) к безразмерным величинам (5) и сохраняя только линейные слагаемые, получим основное уравнение

В случае открытых концов трубы давление, в соответствии с (8), должно удовлетворять условиям р(0, £) = р(1, £ = 0). Для задач горения конкретный вид Q( х, £) зависит от выбранной модели подвода теплоты, т.е. зависит от искомых параметров задачи.

3. Модель горения

Наиболее сложно смоделировать механизмы связей при возбуждении акустических колебаний горением. В работе [1] в самом общем виде даны

(3а)

(6)

Искомые величины выражаются через потенциал скорости ф и Q

и = —

дф

дх ’

(7)

(8)

(9)

(10)

основные принципы построения их. Анализ механизма обратных связей для задач горения капель в акустическом потоке газа дается в работе [2].

Различают подвод теплоты в точках, неподвижных относительно стенок трубы, и подвод теплоты на подвижной поверхности. В работе [3] эти два варианта различаются как локальная и конвективная модели подвода теплоты. Конкретный вид Q( х, г) меняется от большого числа режимных и конструктивных параметров, фигурирующих в расчетной модели. Между акустическими колебаниями и колебанием процесса горения стоит целая цепочка связывающих их процессов. На устойчивость процесса горения влияют время запаздывания тепловыделения, качество распыла топлива, полнота сгорания топлива, коэффициент избытка воздуха и т. п. Здесь рассматривается механизм поддержания колебаний столба газа за счет линейной связи между возмущениями теплоты и возникающими в трубе продольными скоростями.

Возмущенное значение теплоты представляется в виде

Q = Q(x)q( х, г), (11)

где б(х ) дает распределение подведенной теплоты по длине трубы, а ^(х, г) определяет амплитуду возмущения.

Связь между q и и в общем случае представляется в виде

А ()= ¥Ь2 (и), (12)

где Ь,1, ^2 - линейные дифференциальные операторы, зависящие от времени. Величина У в данной работе далее полагается постоянной.

3.1. Локальная модель. Рассмотрим три варианта этой модели, часто используемые при теоретическом анализе [3]. Гармоническое решение системы, не зависящее от начальных условий, отыскивается в виде

q(x, г)= q(x) е1'®*, и(х, г)= и(х) е1'®*. (13)

Линейная связь между этими возмущениями представляется в виде

q = У1 (о) и, I (о) = Iл (о)+ /// (®), (14)

где 1 (о) - комплексное число. Значения Iл (о) и II (о) определяют выбранную связь между параметрами в первом уравнении (14). В работе [3] она трактуется как локальная модель.

Конкретный вид q, Iл и II для трех упомянутых вариантов в случае

локальной модели представляется в таблице 1. Такое задание полностью определяет соотношение (12) для выбранной модели.

Таблица 1

Связь между q и и в случае локальной модели

Вариант = IR II

1 q = Уи( - т) С08(ют) - 8ш(<Вт)

Вариант = q 1R 11

2 z(dq/dt)+ q = Yu 1 /[1 + (ют)2 ] _ ют /[1 + (ют)2 ]

3 q = Y [u + T(dq/dt)] 1 ют

3.2. Конвективная модель. Более сложное соотношение между теплотой и скоростью задается в конвективной модели горения. В этом случае теплота удовлетворяет уравнению

Л Л

+ ис(х,t)-^ = 0, q(xs ,t)= qs , (15)

dt ex

где uc - заданная конвективная скорость, она отличается от скорости частиц газа; xs - координата источника теплоподвода. В соответствии с этой моделью все параметры вниз по течению от сечения х > xs определяются конвективной

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

Текущее время 0, которому соответствует координата xs, находится из решения уравнения dt = dx/uc в предположении uc = const:

0 = t - (x - xs )/uc .

Решение уравнения (15) с учетом значения 0 имеет вид

/ x _ x Л

q(x, t)= qs t----— H(x _ xs ), (16)

I uc

в котором H - единичная ступенчатая функция. Из представления (16) видно,

что поступление теплоты начинается с момента времени t >0, т.е. когда

координата x> xs.

Полагаем, что между qs и us задана линейная дифференциальная связь вида (12)

L1 (qs )= Y(t )L2 (us ). (17)

Периодическое решение отыскивается в виде

qs(t )=qseKOt, us(t )=usel(at. (18)

Тогда соотношение между возмущениями теплоты и скорости

представляется в виде

„ г/ \ -гю(х-х_)/1

Че = УІ(©) ыяв ' •>

Уравнение (19) имеет тот же вид, что и (14), но только со степенным множителем, который дает сдвиг по фазе между возмущениями. Выбор расчетного варианта зависит от конструктивных особенностей камеры и самого процесса горения. Для конвективной модели три варианта связи (17) между возмущениями теплоты и скорости приведены в таблице 2.

Таблица 2

Связь между qs и в случае конвективной модели

вариант N ъ ІЯ ІІ

1 и 1 вч £ N ъ С0в(©т) - 8Іп(©т)

2 )+= Уые 1/ 1 + (©т)2 -©т/ 1 + (©т)2

3 Ч* = У [ы* +т(д*М)] 1 ©т

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

4. Метод решения

Решение уравнения (6) отыскивается в виде суммы

<» Я

♦(хІ) = £ /т () Vт (х); І) = 2 8 т (^'т (х)

т=0 т=0

или (20)

№

Че (-^ І)= 2 8т С)ч'т (х* ).

т=0

Функции V п (х) являются решением однородного уравнения

УП +0 п V п = ^

*2

(21)

где Оп - квадрат собственной частоты, связанный с акустической задачей для трубы, а уп (х) являются собственными функциями, удовлетворяющими условиям ортонормировки. В случае открытых концов трубы функция имеет вид уп(х)= 428ш(Опх) и Оп = п п (п = 1, 2, 3, ...). Коэффициент 42 вводится ввиду выполнения условия нормировки.

В результате специального выбора базисных функций у п (х) система уравнений для определения амплитуд /п (*) и gn (*) находится методом Галеркина. После подстановки (20) в уравнения (6) - (10) и интегрирования от х = 0 до х = 1 с учетом (21) получим

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

Система уравнений (22) в общем случае является взаимосвязанной. Это видно из представления Qn в виде интеграла (23). Ясно, что для сложной

Для исследования задачи устойчивости используем соотношение (22) в простом случае, когда взаимодействиями между модами колебаний пренебрегаем, т.е. при т = п. В результате такого упрощения можно проанализировать устойчивость решения уравнения (22) по отношению к геометрическим и термодинамическим параметрам задачи. Гармоническое решение отыскивается в виде

где /п и gn - постоянные величины. Величина гап является комплексной частотой, по реальной части которой можно судить об устойчивости решения дифференциального уравнения (22). Подставляя (24) в (22) с учетом принятых допущений, получим

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

При принятых выше ограничениях относительно взаимодействия мод колебаний последнее слагаемое в (25) трактуется как величина второго порядка

т=0

или

(23)

функции

вычисление соответствующих интегралов возможно только

численно.

(24)

(25)

или

(26)

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

(<»п -ОпХ<»п + 0п) = (У- і)впп /(©п). (27)

В первом приближении, пренебрегая слагаемыми в правой части (27), получим © п = ±Оп. Второе приближение для комплексной частоты

© п =± Эп + ^ п (28)

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

Эп =О п + (^О ^ \ря,пп/я (о п Х В/,пп// (о п )], (29)

20 п

Сп =^^01Х^ Р,пп// (Ои )-»/,пп/я (Ои )].

2П п

5. Обсуждение результатов

Распределенная теплота по длине трубы задается в форме трапеции

6 (х )= 60 Н (X-А^ - Н (X-А 2)] + н (X-А 2) -

[Л2 - А1

(30)

-Н(х -- з) + Т-4—— р (х -- з) - Н(х -- 4 )]1,

Л4--3 ]

“3) + -4-Х р( х--з) - Н( х--4)], (31)

4 - -3

где 60, А1, А2, А3,А4 — заданные константы, причем А 1< А2< Аз< А4. Зона

подвода теплоты располагается в интервале А1 < х < А4. Придавая разные

значения А; в (31), можно получить распределение теплоты по длине трубы в

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

Выражение для коэффициента ^1 в случае локальной модели найдем в результате подстановки (31) в первое уравнение (26). Опуская возникающие в ходе вычислений интеграла (26) громоздкие выкладки, запишем окончательную формулу:

51 = 4Й7 воГ' ™("1 1 В 2 +Я1 ) )] +

+ 2$т[° 1 (я 3 +Я 2 )]8т[° 1 (Я 3 -Я 2 )]-^[р 1 (я 4 +Я 3 )] ' О ЯЯ 4 Я—). (32)

О 1 Iя 4 -Я 3 /

Аналогичная формула для коэффициента затухания в случае конвективной модели находится в результате подстановки (31) во второе уравнение (26):

Цп.

Ап

1

81П

1

+ -

В2п

1

В

2п

В

+

1п 1

В1п

1

В

2п

В 1п

(

э1п £ 0

V

(

С0Э £ 0

V

(

С0в £ 0

V

(

э1п £ 0

V

(

э1п £ 0

V

1п

Я2 + Я1 1 81п(В1и (Я2 -Я1) 1 2)

В1пВ2 -Я 1)/ 2

2п

В2п

1п

1п

2п

Я 2 +Я1

2

81п(В2и (Я2 -Я1)|2)

В2п (Я2 -Я 1 )| 2

+

Я3 +Я2 >

2

Я 3 + Я 2

Э1П

0п + В

2п

Я3 - Я2

2

Я4 + Я3 2

Я 4 + Я 3 2

Э1П

0п + В

1п

Я3 - Я2

+

в1п(В1п (Я4 -Я3 )|2)

ВЫ (Я4 -Я3)/ 2 81п(в2и (Я4 -Я3 ) 1 2)

В2п (Я4 -Я3 )|2 Здесь для краткости записи введены обозначения

У-1 „ /_ \ „ _ л, О п (1 - ис)

- бо^и С08(°пх.,); 0п = °п - Оп ; В1п =

Ап

В2п =

О п (1 + ис )

(33)

(34)

Выражения (32) и (33) определяют величину и знак коэффициента затухания в зависимости от геометрических параметров задачи и от выбранного закона подвода теплоты (32). Из соотношения (29) можно найти поправку к частоте колебаний столба газа.

5.1. Пример расчета по локальной модели. На рисунках 1-4 представлены результаты расчетов по формуле (32). Расчеты проведены для варианта 1 из таблицы 1. Если в (32) положить Я 2 = Я 3, то распределение теплоты вдоль длины трубы задается в виде треугольника. В работе [3] для этого случая при Я1 = 0 построены значения коэффициента затухания в зависимости от Я 2 при различных значениях Я3 . Если разница Я3 - Я2 ^ 0,01, то приведенные кривые на рис.1 практически совпадают с аналогичными зависимостями в работе [3].

2

2

1

2

с

с

и

с

Исследуем изменение коэффициента £ в зависимости от ширины верхнего среза импульса X3 -X2, сохраняя параметры X1 и X4 постоянными. Для X 4 = 0,9 и при X з — X 2 = 0,01 (рис. 1) условие устойчивости выполняется при А2 ^ 0,55. С увеличением разницы X3 — X2 до 0,1; 0,2 и 0,5 (последовательно рисунки 2,3 и 4) значения X2, которым соответствует устойчивое поведение газа (£ 1 > 0), смещаются влево и принимают значения X2 = 0,48; 0,42 и 0,27, соответственно.

Рис. 1. Коэффициент затухания колебаний для Рис. 2. Коэффициент затухания колебаний

локальной модели для локальной модели

при X1 = 0,0 и X з — X 2= 0,01 при X 1=0.0 и X 3 — X 2 = 0,1

Формы зависимости теплоты 2 от координаты в рассматриваемом случае для X 4 = 0,9 представлены на рисунках 5-8. Из анализа этих рисунков видно, что увеличение площади между ломаной линией и осью 0х расширяет границы устойчивости колебаний столба газа по координате X2 .

С уменьшением максимального значения X 4 до 0,8 уменьшается и соответствующая площадь под кривой теплоты (рис. 9-12).

Рис. 3. Коэффициент затухания колебаний для Рис. 4. Коэффициент затухания колебаний

локальной модели при X1 = 0,0 и X 3 — X 2 = 0,2

для локальной модели при X1 = 0,0 и X 3 — X 2 = 0,5

1 п О

о

Рис. 5. X2 = 0,56, X 3 = 0,57, X 4 = 0,9 Рис. 6. X 2 = 0,48, X 3 = 0,58, X 4 = 0,9

1

О

Рис. 7. X 2 = 0,42, X 3 = 0,62, X 4 = 0,9

1

О

Рис. 9. X 2 = 0,61, X 3 = 0,62, X 4 = 0,9 © Проблемы энергетики, 2004, № 3-4

Рис. 8. X 2 = 0,26, X3 = 0,76, X 4 =0,9

1

О

Рис. 10. X 2 = 0,6, X 3 = 0,6, X 4 = 0,8

о

о

о

о

Для Аз -А2 = 0,5 устойчивый вариант процесса поддерживается при X2 ^ 0,28. Правая граница теплоподвода в этом варианте становится практически вертикальной (А3 » А4) (рис. 12).

Рис. 11. А 2 = 0,44, А 3 = 0,64, А 4 = 0,8 Рис. 12. А 2 = 0,28, А 3 = 0,78, А 4 = 0,8

5.2. Пример расчета по конвективной модели. Рассмотрим частный случай представления (31), в котором в фигурных скобках оставлены второе и третье слагаемые. С учетом переиндексации А; имеем

б (х) = 60 [(х - А1)- Н(х - А2 )].

В решении (33) для этого случая остаются только третье и четвертое слагаемые, в которых значения А 3 и А 2 заменяются на А 2 и А1 соответственно. На рис. 13 представлена зависимость коэффициента затухания в зависимости от конвективной скорости. В соответствии с таблицей 2 расчеты проводились для варианта 1. Из рисунка видно, что с уменьшением конвективной скорости (ис < 0,18) величина £ начинает часто менять знак. Для А1 = 0,4, А 2 = 1,0 устойчивые решения находятся в пределах 0,19 < ис < 0,84; для А1 = 0,5, А2 = 0,90 - в пределах 0,31 < ис < 0,94 . Если А1 = 0,60, А2 = 0,80, то интервал устойчивых решений находится в пределах 0,33 < ис < 0,97. С уменьшением разницы (А2- А1) при фиксированном значении суммы (А2 +А1 )= 1,4 левая граница интервала смещается по оси 0х вправо, при этом длина его незначительно уменьшается.

На рис. 14 построены зависимости £ от А = (А2 +А1)/2. Все кривые построены с учетом выполнения условия 0 < 0 < А1 < А2 < 1,0, поэтому разные кривые имеют свое начало и конец. Из рисунка видно, что устойчивые колебания наблюдаются при А>0,64. С увеличением А2 -А 1 = 0,1 положительные значения £ находятся в интервале 0,64 < х < 0,94. Для А 2 - А 1 = 0,2 интервал устойчивых решений имеет границы 0,65 < х < 0,89. Минимальный интервал устойчивых решений наблюдается для А 2 - А 1 = 0,4 . При А 2 - А 1 > 0,4 решения неустойчивы.

Рис. 13. Зависимость коэффициента затухания Рис. 14. Зависимость коэффициента

от скорости uc для конвективной модели затухания от X для конвективной модели

6. Выводы

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

Summary

A mathematical model of unsteady velocity sensitive thermoacoustic interactions in a heated duct (generalized Rijke tube) is developed. The model is constructed in a general way which is capable of employing any linear velocity sensitive thermoacoustic response model. A solution methodology is then constructed based on modal analysis. A number of rather general closed form solutions are presented and some of the results are used to qualitatively explain the behaviors of Rijke tubes.

Литература

1. Раушенбах Б.В. Вибрационное горение. - М.: ГИФМЛ, 1961. - 500 с.

2. Натанзон М.С. Неустойчивость горения. - М.: Машиностроение, 1986. -248 с.

3. Hyun-Gull Yoon, John Peddieson Jr., Kenneth R. Purdy. Mathematical modeling of a generalized Rijke tube // International Journal of Engineering Science. -1998. - № 36. - рр. 1234-1264.

4. Юдаев Б.Н. Теплопередача. - М.: Высшая школа, 1981. - 320 с.

Поступила 22.09.2003