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

УДК 621.436.052

Ю. А. Гришин, Н. В. Х а з о в

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

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

Ключевые слова: поршневой двигатель, газовоздушный тракт, газодинамика, метод характеристик, моделирование течений.

В последнее время при создании систем газообмена поршневых двигателей широкое применение нашли компьютерные программы моделирования газодинамических течений, так называемые CFD-комплексы. Среди них выделяется большой класс комплексов, выполняющих расчеты нестационарных процессов и конструкторскую оптимизацию газовоздушных трактов двигателей в одномерной постановке. Таковыми являются Boost (AVL), Wave (Ricardo Software) и др. В их основе лежат те или иные численные методы газовой динамики.

Как один из наиболее эффективных методов моделирования одномерных течений выделяется метод характеристик с плавающей сеткой (МХПС), позволяющий получать точное фазовое и амплитудное соответствие расчетных процессов действительным [1-4]. Одним из недостатков данного метода, ограничивающим его применение [5], является отсутствие в его алгоритме интегральных балансовых соотношений. Это приводит к неточностям расчета таких ключевых параметров, как объем воздуха, попавшего в цилиндры двигателя, и располагаемая энергия газов перед турбиной.

Новая версия МХПС, представленная в работах [6, 7], позволяет обеспечить выполнение законов сохранения, т.е. является консервативной. В настоящей работе приведено краткое описание метода и результаты тестовых расчетов, соответствующих экспериментальным данным, полученным на специальной установке.

Рассмотрим построение расчетной сетки и определение газодинамических параметров на примере ячейки d, изображенной на рис. 1. В ячейках a,b и других ячейках предыдущего ряда значения газодинамических параметров уже известны так же, как и координаты узлов М, С, N и т.д. В настоящей работе приведен уточненный, квазинеявный способ определения угловых коэффициентов (скоростей фронтов) wac и wbc отрезков характеристик AC и BC. Для этого, используя

Рис. 1. К построению ячейки d

известные соотношения распада разрыва [3] или газодинамические функции нестационарного течения [7, 8], находят предварительные значения скорости u и скорости звука a в зонах da и db:

Щ + Ща aa + ada

wac = (u - a)Ac =

2 2 + + (1)

Щ + Щъ +

wвc = (и + а)вс =-2--1--2-'

Аналогичным образом определяют угловые коэффициенты отрезков МА и БЫ, что позволяет в итоге получить координаты гА, хА и гв, хв. Запишем значения массы Ма, количества движения Ка и полной энергии Еа, которые будут переходить из ячейки а в ячейку d через фронт АС при единичной площади проходного сечения:

Ма = ра (иа - W АС) & А - ^С) ; (2)

К а = Ма Ща + Ра (рА - ^ );(3)

Еа = Ма ва + Ра Ща (рА - ^), (4)

2 2 _ и р и2

здесь р, р, в = вуТ +--= —----1---плотность, давление и

2 (к - 1)р 2

удельная энергия. Подобным образом рассчитывают значения Мъ, Къ

и Еъ для перехода из ячейки Ь в ячейку d:

Мъ = ръ (иъ + wвc) (рв - гс); (5)

Къ = Мъ иъ - Ръ (гв - гс );(б)

Еъ = Мъ въ - Ръ иъ (г в - гс). (7)

Всего в ячейку d поступит М = Ма + Мъ, К = Ка + Къ и Е = Еа + + Еъ. Сверху ячейку d будут замыкать отрезки характеристик Б А и Б Б. Полученные из нижнего ряда суммарные значения М, К и Е должны быть полностью распределены по ячейкам / и д следующего ряда с помощью соотношений (2)-(7). Тем самым в контуре САББ

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

При этом для массы можно записать

М = Мщ + МЛд = ра(хол - хл) + Ра(хв - хив) =

= Ра [(хв - хл) + (хил - хив)]. (8)

Для балансовых соотношений количества движения и энергии получаем соответственно

К = Мщиа - ра^и - 1ил) +

+ МЛдиа + ра^и - 1ив) = Млиа - ра^ив - ^ил); (9)

E = M,

ff

Pd

u

+ — (k - l)pd 2 _

- PdUd(tD - tDA) + Mf

dg

Pd

u

(k - l)pd + 2 _

+ PdUd(tD - tDB) =

= Md.

Pd

u

21

(k - l)pd + 2 .

- PdUd(tDB - tDA). (10)

Обозначая хв - хл = Ах, 1в - Ьл = Ьив - Ьил = Аt и учитывая, что БлБв — это отрезок линии тока с угловым коэффициентом в виде искомой скорости иа = (хил-хив)/(^ил-ив) = (хил-хив)/(-А£), уравнения (8)-(10) можно переписать в виде:

М = ра(Ах - иА); (11)

К = Миа - рА; (12)

E=M

Pd

ud + f

- pdUd At,

(13)

(k - l)pd

что позволяет определить pd, ud и pd.

При построении сетки может получиться At = 0. Тогда из уравнения (11) сразу определяем pd, из формулы (12) — ud, затем из выражения (13) — pd. Для общего случая при At = 0 с использованием обозначений i = K/M, е = Е/М, 8 = Ах/At из системы уравнений (11)—(13) получаем формулу для определения скорости ud:

8 + ki Ud = ~k+l ±

'8 + ki

2

(k - l) е + 8i k+l .

(14)

Теперь из уравнений (11) и (13) легко найти ра и ра, а затем и скорость звука аа = \/кра/ра, необходимую для определения скоростей фронтов БЛ и ББ (см. (1)). При расчете иа для выбора знака перед радикалом в уравнении (14) необходимо использовать единственное контрольное значение скорости, получаемое из уравнения (12) при

2

АЬ = 0, т.е. ud = 1. Из двух значений иа по формуле (14) нужно выбрать одно, наиболее близкое к значению 1.

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

Данная схема была реализована на ЭВМ и проверена на выполнение балансов по массе, количеству движения и энергии потоков, поступивших и вышедших из трубопровода за произвольный промежуток времени. Тем самым была подтверждена консервативность данной версии метода характеристик с плавающей сеткой.

На следующем этапе проверки на адекватность данной схемы было выполнено численное моделирование и проведено сравнение с реальным нестационарным процессом. Известно, что для проверки численных моделей наиболее точные результаты можно получить при исследовании формирования и движения уединенных волн конечной амплитуды. Часто при этом применяют ударные трубы. Однако известно, что при диагностике крутых ударных фронтов трудно обеспечить точность записи давления за фронтом из-за физических особенностей работы датчиков [9]. Наиболее точную запись можно получить при гладкой, скругленной форме импульсов. Кроме того, важно отметить, что именно такую форму имеют импульсы, поступающие из органов газообмена в каналы впускных и выпускных систем реальных поршневых двигателей. Необходимые эксперименты были проведены на кафедре "Поршневые двигатели" МГТУ им. Н.Э. Баумана с помощью специального генератора уединенных волн (рис. 2).

Исходный газодинамический импульс скругленной формы (уединенная волна) получался при плавном открытии и последующем закрытии клапана, соединяющего цилиндр высокого давления воздуха с длинной гладкостенной трубой. Импульс записывался датчиком давления р!, установленным на длине трубы Ь! = 0,2 м, где обеспечено присоединение потока к стенкам канала после отрыва за клапаном. Для записи импульса, отраженного от открытого конца трубы, целесообразно использовать датчик р2 (Ь2 ~ 0), где амплитуда в результате взаимодействия с закрытым концом увеличивается (клапан в это время уже закрыт).

При численном моделировании реальных процессов необходимо учитывать трение и теплообмен. Здесь можно использовать представления Дарси-Вейсбаха и Ньютона — так называемый гидравлический

Рис.2. Cхема экспериментальной установки для исследования уединенных волн

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

\_pudx — (- + pu2) dt] = —

pedx — pv ( e + — ) dt p

dxdt ;

Xpu2 ~2D

4a (T — Tk) D

(15)

dxdt, (16)

где Л и а — коэффициенты трения и теплоотдачи; Б и Тк — внутренний диаметр и температура стенки трубы. С использованием конечно-разностной записи правых частей уравнений (15) и (16) значения К и Е, подставляемые в формулы (12) и (13), следует уменьшить соответственно на величины

Xpaua

2D

(Хс — XAv )(tA — tc) +

Л

Л^к

2D

(xBv — XC )(tB — tC) =

= 2D \mau2a(tA — tc)2 + mb^(tB — tc)2] ; (17)

4a[-a/(k — l)pa — Tk ]

D

(xc — xav )(tA — tc) +

+

4a[pb/(k — l)pb — Tk ]

4a

D

D

-a

(k — l)pa

-b

Tk

(xBv — xc )(tB — tc) = (ua — Wac )(tA — tc )2 +

+

(k — l)pb

Tk

(wbc — ub)(tB — tc)2

(18)

Очевидно, что при этом необходимо знание коэффициентов Л, а. Наиболее важен учет трения, поскольку оно наблюдается даже при течении ненагретого газа. Тонким моментом при моделировании результатов описываемого эксперимента являлось задание значения коэффициента Л, который, как известно, зависит от числа Рейнольдса. На первой стадии моделирования его определяли для каждого расчетного шага по известной формуле Блазиуса:

0,316

Л=

Re'

0,25 •

(19)

На рис. 3 показаны соответствующие результаты расчета МХПС отраженной волны с учетом трения при запуске 30 характеристик, т.е. при 30 расчетных ячейках по длине трубы, и результаты экспериментального осциллографирования.

Как видно из рис.3, в отраженном импульсе имеется 5 %-ная погрешность моделирования по амплитуде давления. Это подтверждает точку зрения о неточности применения стационарной зависимости для нестационарных процессов. Следует отметить, что такая ошибка

Рис.3. Результаты осциллографирования волнового процесса (сплошная линия) и расчетные результаты (точки), полученные МХПС c использованием фомулы Блазиуса (19)

при определении давления приводит к гораздо более значительным ошибкам при определении расходов — до 15 %, что недопустимо при расчетах газообмена в поршневых двигателях.

Для коэффициента трения при нестабилизированном течении на начальных участках труб известна также формула [10]

0,344

Л = (Rex/D )0,2, (20)

где х — длина участка стабилизации течения.

Однако подстановка этого выражения в программу расчета привела лишь к незначительному улучшению результатов по сравнению с предыдущими. Поэтому было принято решение о проведении специальных исследований по уточнению коэффициента Л. Для этого с помощью описанного ранее генератора были проведены многочисленные эксперименты на длинных трубах (L = 3,5 м) с установкой второго датчика давления на длине L2 = 1,55м (рис.4). Эту длину определяли специальным образом. С одной стороны, она должна быть как можно больше, чтобы повысить точность оценки снижения амплитуды импульса за счет трения, с другой — фронты импульса не должны успеть превратиться в ударные, что привело бы к снижению точности записи давления. Таким образом, влияние трения на снижение амплитуды оценивалось на длине L2 — Li = 1,35 м. Удлинение труб до 3,5 м обеспечивало "чистую" запись импульсов вторым датчиком, так как гарантировалось отсутствие волн, отраженных от открытого конца.

Внутренний диаметр гладких труб из коррозионно-стойкой стали варьировался от 20 до 42 мм, испускались как положительные (а) так и отрицательные (б) импульсы давления различной амплитуды в диапазонах соответственно от 0 до +80 кПа изб. и от 0 до —50 кПа изб. Такое ограничение диапазонов амплитуд обосновано, во-первых, тем обстоятельством, что в системах газообмена поршневых двигателей

PI 1550 PI

200

П

3500

P

P

а

6

Рис. 4. Схема эксперимента для определения степени затухания исходных импульсов вследствие трения

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

Результаты экспериментов показаны на рис.5. Из рисунка видно, что амплитуда импульсов на мерном участке трубопровода заметно уменьшается как в области положительных, так и в области отрицательных импульсов.

Затем были проведены расчеты этих импульсов МХПС с подбором коэффициентов трения Л, которые обеспечивали снижение амплитуды р2 до экспериментальных значений. На рис. 6 точками показаны обработанные результаты в системе координат ^ Яе — ^(100Л), т.е. была решена обратная задача — получены значения Л при различных числах Рейнольдса (см. рис. 6). Для обработки были использованы также результаты из двух публикаций, где расчеты импульсов проводились неконсервативными версиями МХПС: по известной схеме Массо [3] и с использованием более совершенной схемы без третьего характеристического семейства [4].

Из рис. 6 видно, что практически для всей исследуемой области нестационарного течения по числу Яе значения Л намного превышают соответствующие значения, которые имеют место при стационарном течении в трубах. Во многих работах, в частности в работах [9, 11, 12], отмечается более сложный характер вязких эффектов при нестационарном течении, которые приводят к затуханию волновых процессов. Помимо увеличения напряжения трения на стенках труб вводится предположение о возникновении на стенках слабых поперечных волн, усиливающих эффект диссипации.

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

Л - Re12 •

31600

(21)

Рис. 5. Влияние трения на амплитуды идущих импульсов

lg(100X)

1,2

1,0

0,8 0,6 0,4 0,2

у 1 1мпул 1 [ЬСЫ

\ □

\ \ , к f

V/ § \ X

Стационарное течение 1 N. [

| | [

2,5 3,0 3,5 4,0 4,5 5,0 ^(Яе)

Рис. 6. Зависимости коэффициента трения Л от числа Рейнольдса Re при стационарном и нестационарном течениях. Расчеты импульсов МХПС:

по схеме консервативной версии данной работы х — импульс сжатия, о — импульс разрежения; по схеме работы

[4] +--импульс сжатия, □ — мпульс

разрежения; по схеме Массо [3] — А

Подставляя выражение (21) в расчетную программу, получаем практически полное совпадение и по фазе, и по амплитуде результатов эксперимента и моделирования отраженного импульса при запуске даже малого числа характеристик (порядка 30), что изображено на рис. 7.

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

Рис. 7. Результаты осциллографирования волнового процесса (сплошная) и расчетные результаты, полученные c использованием формулы (21): новая версия МХПС (точки), метод РПР (штриховая кривая)

Для сравнения на рис. 7 показаны результаты расчета, выполненного с учетом трения по зависимости (21) известным консервативным численным методом распада разрыва (РПР) [2, 3] с тем же числом ячеек. Эти расчеты показали значительное искажение формы импульса, особенно его заднего фронта; по амплитуде неточность составила 20%. В целях получения результата с нужной точностью оказалось необходимым увеличивать число расчетных ячеек до 600, при этом время вычислений возрастает примерно в 10 раз по сравнению с изложенной новой версией МХПС.

По результатам данной работы можно сделать следующие выводы.

1. Представлены основные расчетные соотношения новой версии численного метода характеристик с плавающей сеткой, который обеспечивает выполнение интегральных законов сохранения.

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

3. Показана неточность, которая может быть получена при расчетах импульсов в случае использования известных соотношений для коэффициента трения Л.

4. С помощью специальных экспериментальных исследований определена степень затухания импульсов в широком диапазоне амплитуд и чисел Re.

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

6. Численные расчеты прошедших и отраженных импульсов, выполненные с использованием полученной зависимости Л = Л(Re), показали хорошее согласование с результатами соответствующих экспериментов.

7. Показаны преимущества разработанной консервативной версии метода характеристик с плавающей сеткой по сравнению с методом распада разрыва, в частности по времени счета.

СПИСОК ЛИТЕРАТУРЫ

1. Benson R., Galloway S. An experimental and analytical investigation of the gas exchange process in a multicylinder pressure-charged two-stroke engine // Proc. In. Mech. Eng. - 1968. - V. 183. Pt. 1. No. 14. - P. 253-267.

2. ПирумовУ. Г., Росляков Г. С. Численные методы газовой динамики. -М.: Высш. шк., 1987. - 323 с.

3. Р у д о й Б. П. Прикладная нестационарная гидрогазодинамика: Учеб. пособие. -Уфа: УАИ, 1988. - 184 с.

4. Гришин Ю. А. Версия метода характеристик с плавающей сеткой // Математическое моделирование. - 2002. - Т. 15. - № 8. - С. 3-8.

5. РetersB., GosmanA. Numerical simulation of unsteady flow in engine intake manifolds // SAE Paper № 930609.

6. Гришин Ю. А. Консервативный метод характеристик с плавающей сеткой // Тез. докл. V Междунар. конф. по неравновесным процессам в соплах и струях (NPNJ), Самара, 5-10 июля 2004 г. - М.: Вузовская книга, 2004. - С. 79-80.

7. Гришин Ю. А. Численное решение практических задач газовой динамики в поршневых двигателях // Изв. ТулГУ Сер. Автомобильный транспорт. - 2005. -Вып. 9. - С. 173-179.

8. Гришин Ю. А., К р у г л о в М. Г., Р у д о й Б. П. Газодинамические функции для расчета нестационарных течений газа // Изв.вузов. Машиностроение. -1977. -№ 3. - С. 52-56.

9. Ударные трубы: Сб. стат. / Под ред. Х.А. Рахматуллина и С.С. Семенова. -М.: Изд-во иностр. лит-ры, 1962. - 700 с.

10. И д е л ь ч и к И. Е. Справочник по гидравлическим сопротивлениям / Под ред. М.О. Штейнберга. - М.: Машиностроение, 1992. - 672 с.

11. Гришин Ю. А., Григорьев М. А., Пудовкин И. Ю.,ХазовН. В. Исследование влияния трения при расчетах нестационарного течения в протяженных каналах // Тез. докл. V Междунар. конф. по неравновесным процессам в соплах и струях (NPNJ), Самара, 5-10 июля 2004 г. - М.: Вузовская книга, 2004. - С. 81-82.

12. Основные результаты экспериментов на ударных трубах / Под ред. А. Ферри. - М.: Гос. изд-во литературы по атомной науке и технике, 1963. -442 с.

Статья поступила в редакцию 10.07.2008

Юрий Аркадьевич Гришин родился в 1947 г., окончил в 1971 г. Уфимский авиационный институт им. С. Орджоникидзе. Д-р техн. наук, профессор кафедры "Поршневые двигатели" МГТУ им. Н.Э. Баумана. Автор более 160 научных работ в области газовой динамики и расчетно-экспериментальных исследований двигателей.

Yu.A. Grishin (b. 1947) graduated from the Ufa Aviation Institute n.a. S. Ordzhonikidze in 1971. D. Sc. (Eng.), professor of "Reciprocating Engines" department of the Bauman Moscow State Technical University. Author of more than 160 publications in the field of gas dynamics and computational and experimental studies of engines.

Николай Владимирович Хазов родился в 1982 г., окончил в 2005 г. МГТУ им. Н.Э. Баумана. Аспирант кафедры "Поршневые двигатели" МГТУ им. Н.Э. Баумана. Автор 2 научных работ по численным исследованиям в области газовой динамики двигателей.

N.V. Khazov (b. 1982) graduated from the Bauman Moscow State Technical University in 2005. Post-graduate of "Reciprocating Engines" department of the Bauman Moscow State Technical University. Author of 2 publications in the field of numerical studies in the field of gas dynamics of engines.