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

Математические проблемы расчета многозонных моделей рабочего цикла

ДВС с искровым зажиганием

к.т.н. доц. Апелинский Д.В., доц. Шендеровский И.М., к.т.н. доц. Руновский К.С.

Университет машиностроения 8(495) 223-05-23 (доб. 1467), ареИтИутатЩтаИги

Аннотация. В статье рассмотрены особенности математических моделей рабочих процессов в ДВС, не позволяющие выполнить приведение обыкновенных дифференциальных уравнений (ОДУ) этих моделей к форме Коши. Предложен абсолютно устойчивый алгоритм численного интегрирования систем ОДУ, не разрешенных относительно производных, использующий сплайны Эрмита и методы нелинейного программирования.

Ключевые слова: метод контрольных объемов, многозонная модель, модель Близарда - Кекка, уравнения баланса массы и энергии, баротропное течение, задача Коши, системы обыкновенных дифференциальных уравнений, не разрешенных относительно производных, пространство джетов, сплайны Эрмита, среднеквадратичная погрешность, нелинейное программирование. Для того чтобы рационализировать процедуры конструирования новых экологически чистых отечественных ДВС и средств электронного управления ими, кафедрой «Автомобильные и тракторные двигатели» Университета машиностроения проводится разработка системы математических моделей и программного обеспечения для компьютерного моделирования рабочего цикла этих двигателей. Методической основой этой работы являются т.н. детерминистические модели физико-химических процессов, протекающих в цилиндрах и газовоздушных трактах ДВС, т.е. модели, использующие классические законы газовой динамики, термодинамики и химической кинетики, сформулированные на языке дифференциальных и интегральных уравнений и описывающие эволюцию рассматриваемой системы и ее частей.

Существуют два метода конструирования таких моделей.

В первом из них, возникшем в 60-х годах XX века на основе работ проф. И.И. Вибе и Г. Вошни по моделированию внутрицилиндровых процессов, рабочее тело представляется как совокупность термодинамически однородных зон, рассматриваемых как открытые термодинамические системы. Для каждой из зон составляются дифференциальные уравнения балансов массы и энергии, а также уравнение состояния рабочего тела и уравнения, описывающие изменение объема зоны и площади ее поверхности. Недостатком этого метода является отсутствие адекватного описания газодинамических процессов в цилиндрах и системе газообмена двигателей.

Во втором методе, возникшем на базе численных методов газовой динамики и теплофизики в 80-х годах XX в., рабочее тело в цилиндре и газовых трактах двигателя рассматривается как совокупность т.н. контрольных объемов (КО). Для каждого КО составляются уравнения баланса массы, энергии, импульса и других физических и химических экстенсивных параметров. Применение указанной процедуры к бесконечно малым КО позволяет получить систему дифференциальных уравнений в частных производных, описывающих движение многофазной химически реагирующей материальной среды. К сожалению, вычислительная сложность процедуры решения этих уравнений препятствует их использованию в качестве интеллектуального инструментария инженера при проектировании и модернизации ДВС. Для сокращения объема вычислительной работы переходят к конечным контрольным объемам (ККО), для которых вышеуказанные уравнения баланса записываются в виде интегралов по выделенным объемам конечной величины и их поверхностям. Введение аналитической аппроксимации рассматриваемых функций внутри и на поверхности ККО позволяет свести задачу моделирования к решению системы обыкновенных дифференциальных уравнений (ОДУ), что требует значительно меньших компьютерных ресурсов.

Объединив достоинства обоих методов, в теории ДВС пришли к концепции многозон-

Серия 1. Наземные транспортные средства, энергетические установки и двигатели.

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

Например, описывая процесс сгорания в бензиновых двигателях по методу Близарда-Кекка [1], в качестве ККО выбирают объем надпоршневого пространства, характеризуемый давлением p, объемом Vcyl, массой m, кинетической энергией регулярного течения K и

кинетической энергией турбулентности к . Рабочее тело подразделяют на 3 зоны:

• зону несгоревшей (unburned) части заряда, характеризуемой массой mu, объемом Vu, температурой Tu, удельными плотностями внутренней энергии Uu и энтальпии hu, удельной теплоемкостью cp u = dhujdT и газовой постоянной Ru ;

• зону сгоревшей (burned) части заряда, характеризуемой массой mb, объемом Vb, температурой Tb, удельными плотностями внутренней энергии ub и энтальпии h, удельной

теплоемкостью С,

dhj dT

,р ь — и газовой постоянной яь,

зону вовлечения (еп1хаттеп1:), состоящую из глобул несгоревшей части заряда, вовлеченных турбулентными пульсациями в зону сгоревшей части заряда, и характеризуемой параметрами те, Уе и Те = Ти,

Модель описывается следующей системой ОДУ: - уравнение баланса тепловой энергии в зоне несгоревшей части заряда

d (muuu) dV dQ

V u u ) p u ^u , w

dt

dt

dt

h

dm

dt

0

уравнение состояния несгоревшей части заряда

1 dp 1 dV 1 dm 1 dT

j u u u

p dt Vu dt

= 0,

mu dt Tu dt

- уравнение баланса тепловой энергии продуктов сгорания

d ( mbub)

+ p

dVb dQb

+

Cb, w -h ^ = 0,

dt dt dt u dt уравнение состояния продуктов сгорания

1 dp 1 dV 1 dm 1 dR,

1 dT

p dt

Vb dt

m

где:

Ru

n

bi'

R dt

dt

Rb dt

T dt

0

-2

mb Z dn^'

dm

v

dt dt

b2X.

(1)

(2)

(3)

(4)

(5)

Я,0 - универсальная газовая постоянная, пь 1 - число киломолей г-го компонента продуктов сгорания, - уравнения кинетики химических реакций в продуктах сгорания dnbl í dmí

m.

dt

+ ("b' - nb! )- VWi?= 0, i = 1.....N

(6)

[е] „ .

здесь: Щ / - число киломолеи г-го компонента в состоянии химического равновесия при заданном давлении р и температуре Ть,

ЖЬ е] - равновесная скорость его образования, вычисляемая по уравнениям химической

кинетики, N - общее число компонентов,

b

- уравнение вовлечения свежей смеси в зону сгорания (1-е уравнение Близарда-Кекка)

dm т Лц

Иш ^ »

^ "^(~Шття + ^ ) = 0. (7)

dt V

здесь: ЛЕ1 - условная площадь фронта пламени, вычисляемая как функция Уь , а также конфигурации камеры сгорания и положения поршня, - среднеквадратическая скорость турбулентности,

- скорость ламинарного пламени в камере сгорания, причем

ЧЬ ЧЬ, 283

Т

( у р

к Ро ;

(4.706ут2 - 4.062ут +1) , (8)

283

SЬ 283 - скорость ламинарного сгорания топливовоздушной смеси при Т = 283°К и Р = Ро,

у т - коэффициент остаточных газов,

ОС и Р - эмпирические константы, зависящие от свойств топлива, - уравнение выгорания свежей смеси (2-е уравнение Близарда - Кеека)

dmь ЧТ , ч

Ь-Чг{те -ть) = 0, (9)

dt X

где: X - микромасштаб Тейлора.

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

- масса т рабочего тела в цилиндре, которая удовлетворяет уравнению:

dm dm dm

—^ + —ь---= 0, (10)

dt dt dt

- среднеквадратическая скорость турбулентных пульсаций: шгт5 = ^2к/(3т) ,

- интегральный масштаб турбулентности: Ь = тт\Ос/ 2, кО^ /41,

- микромасштаб Тейлора: X = ^ 15у Ь/шт , где V - кинематическая вязкость газа,

- уравнение баланса кинетической энергии регулярного течения:

dK dmi„ V2 К dm„

1П 1П

+ Р +--^ = 0, (11)

dt dt 2 т dt

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

Р - С^, (.2)

I V т

в которой I ~ 0.22Ь - линейный размер турбулентности, С - эмпирическая константа, а индексами "т" и "ех" помечены газовые потоки, входящие в камеру сгорания и покидающие

- уравнение баланса кинетической энергии турбулентности:

dk р 2к

Г 7 , / г / Л

су1

V 2

- уравнение объемов:

dt 3¥су1

dm dVc

у су1 dt dt у

3

+ 0.544л/ -кТ = 0, (13)

л' тЬ2

ъ + V _ ^=0. (14)

dt dt dt

Кроме модели процесса сгорания в модель рабочего цикла ДВС должны быть включе-

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

Уравнение баланса тепловой энергии и уравнение состояния при впуске, если не рассматривается процесс испарения топлива в цилиндре, аналогичны соответственно уравнениям (1) и (2), а при выпуске - уравнениям (3), (4) и (5). При рассмотрении процесса испарения топлива на тактах впуска и сжатия, необходимо вместо уравнения (2) использовать уравнение (4). При описании теплообмена между рабочим телом и элементами конструкции цилиндра последние удобно представить в виде совокупности ККО, аппроксимировав температурное поле каждого из них сплайнами Эрмита (см. ниже), и применить к полученной модели стандартное уравнение теплопроводности. В результате имеем систему ОДУ, которая описывает эволюцию температуры и тепловых потоков в узлах сплайнов и решается совместно с уравнениями, описывающими термодинамику рабочего тела. Для описания движения газов в системах впуска и выпуска эти системы моделируются системой, состоящей из М > 1 ККО, адаптированной к особенностям геометрической структуры каналов. Предполагая, что течение баротропно, а уравнение состояния газовой смеси в каждом 1-м ККО описывается уравнением политропы с показателем пг, и, усреднив поле скоростей в поперечных сечениях каналов, получаем следующую рекуррентную систему уравнений:

-Ч а

0, ,

йг

■ + ап

ПгРг

г-1

(Щ 1 - 1)р^

* р , 0

г-1

-п

2 2 , Щ - Щ-1 +

г-1 ( р ^ п,--1

0 1=1

'1

( „ V

Щ -1

-1

К^г-1

= 0,

ц

2п

1 dpi 1 йр^

г ^

йг

йг

ц йР1

+ш Г"

1=2 П1

1 йр 1 йр

1-1

йг

^ -1

йг

(15)

ЩР1

йг

г-1

+ РоП

1=1

V 1 -1 У

Г „ X

w, - w_

г-1

^ г-1 У

= о ,

где: Ц - длина осевой линии г-го ККО, wi - средняя скорость газа в г-м сечении, рг - давление в этом сечении, а0 г - коэффициент количества движения (для турбулентного течения а0г- = 1.07), ар г - коэффициент кинетической энергии, Ъ,р г - коэффициент

гидравлического сопротивления, причем два последних коэффициента приведены к г-му сечению.

Первое из уравнений (15) получено дискретизацией уравнения количества движения, а второе - дискретизацией уравнения неразрывности с учетом политропного уравнения состояния газа в каждом ККО.

При моделировании процесса сгорания в дизеле надпоршневое пространство рассматривается как совокупность ККО, которые конструируются в соответствии с геометрией камеры сгорания. Каждый ККО рассматривается как совокупность зоны несгоревшей части заряда, зоны сгоревшей части заряда, зоны вовлечения и зоны неиспарившегося топлива. Во время впрыска на эту схему накладывается структура топливного факела. Поэтому в уравнение баланса тепловой энергии несгоревшей части заряда добавляются члены, соответствующие энтальпии паров топлива и потери тепла на его испарение. Уравнение состояния несгоревшей части заряда записывается не в форме (2), а в форме (4). Кроме того, вместо уравнения (8), ламинарная скорость сгорания определяется из уравнений химической кинетики процесса с учетом диффузии паров топлива и окислителя в зону фронта пламени. Эти особенности не меняют структуру уравнений, описывающих процесс сгорания, однако значи-

п

г-1

п

тельно увеличивают их количество.

Классические методы численного интегрирования систем ОДУ предусматривают предварительное преобразование их к форме Коши, т.е. к виду:

dxL = fi(t,Xi,x2,...,xN), i = 1,2.....N. (16)

dt

Только после такого преобразования к ним можно применить известные методы численного интегрирования, такие как метод Эйлера, метод Рунге-Кутта, метод Адамса и др. Однако преобразовать уравнения типа (1) - (15) к виду (16) практически невозможно. Даже в тех случаях, когда, упростив коэффициенты некоторых уравнений (например, считая теплоемкости газов постоянными), это преобразование удается выполнить для какой-либо конкретной модели, изменение хотя бы одного уравнения заставляет проделывать всю процедуру преобразования заново и, следовательно, перерабатывать все программное обеспечение задачи моделирования. Кроме того, в процессе решения уравнений (1) - (15) «коэффициенты» при некоторых производных могут стать достаточно малыми, что влечет неустойчивость соответствующего численного метода: численное решение «разбалтывается» из-за накопления погрешностей округления. Такие системы уравнений называются жесткими системами ОДУ. Для работы с жесткими системами ОДУ, возникающими при математическом моделировании электронных цепей, широко используются т.н. неявные методы численного интегрирования ОДУ [2], такие как неявный метод Эйлера, метод Гира - Нордсика и метод Розенброка. Однако программная реализация этих методов жестко привязана к задачам теории цепей, что вызывает значительные трудности при их применении к задачам моделирования рабочего цикла ДВС. К числу недостатков указанных методов следует отнести сложную процедуру выбора длины шага интегрирования и отсутствие механизма фильтрации высокочастотных компонентов решения, которые в рассматриваемой задаче не имеют физического смысла.

Авторами предлагается альтернативный метод численного интегрирования систем ОДУ, не разрешенных относительно производных, свободный от указанных недостатков. Для пояснения логической структуры применяемых алгоритмов рассмотрим ОДУ 1-го порядка, записав его в виде:

f (t, x, X) = 0, (17)

где: x = x(t), а X = dx(t)/dt.

С алгебраической точки зрения уравнение (17) - это уравнение с тремя неизвестными, а его геометрическим образом является поверхность в пространстве переменных |t, x, .x|. В

современной математике набор вида |t, x(t), x(t, состоящий из точки t на оси абсцисс и значений некоторой функции и ее производной в этой точке, называется джетом [3] (струей). Совокупность таких объектов образует пространство джетов /^[t]. Иными словами,

уравнение (17) описывает поверхность в пространстве джетов /^[t], а его решением, удовлетворяющим условию x(0) = x0, является линия на поверхности (17), проходящая через точку пересечения этой поверхности с плоскостями t = 0 и x = x0 .

Естественным способом приближенного описания пространства /^[t ] является применение сплайнов Эрмита. Для этого точками t = tj, j = 0,1,2,...,m,... разбиваем его на систему брусов rij, J = 1,2,____На каждом таком IIj будем аппроксимировать функции x(t)

линейной комбинацией полиномов 3 степени. Ведя вспомогательную переменную I при помощи соотношения t = tJ_1 + hj I, где hj = tj — tj_1, эти полиномы, называемые базисными полиномами Эрмита, можно записать в виде:

2т3 - 3т2 +1, если 0 < х < 1,

ЛДх) ^ —— — — — Л2 (Т) ^

I 0, если т < 0 или т > 1; 10,

т3 - 2т2 + х, если 0 < х < 1,

если х < 0 или х > 1; х3 -х2, если 0 < х < 1,

л , ч I -2т3 - 3х2, если 0 <х< 1, А . . Л3(х) = \ Л4(х) =

I 0, если х < 0 или х > 1; 10, если х < 0 или х > 1.

График функции Л^т)

График функции Л- (ж)

0.2

0.4 0.6 0.8 График функции Л* <с)

Рисунок 1. Поведение базисных сплайнов Эрмита на интервале 0<т<1

При этом для любой функции Х. (г), г -_х < г < г1, записанной в виде:

Х] (г, х10, Х11, х20, Х21 ) = ХюЛ1 (И (г - г}_1 ) + Х1Д Л2 (И (г - г}_1 ) +

+Х20Л3 (И-1 (г - г 1 _1) + х2ХИ } Л 4 (И г1 (г - г}_х),

выполняются соотношения:

йХ(г,...)

Х,-(г,--1>- ■ ■) - Х10,

йг

Х,-(г,->-- ■) — Х20'

йХ.(г,...)

йг

= Х0

(18)

(19)

Поэтому, подставляя в уравнение (17) в качестве искомой функции соотношение (18) и полагая последовательно г = г -_х и г = г1, получаем два алгебраических уравнения с четырьмя неизвестными Х10, Х11, Х20 И Х21. Но по условию задачи для точки г = г0 имеем Х10 = Х(г0 ) = Х0 (начальное условие). Поэтому для формального замыкания задачи на первом

интервале недостает еще одного уравнения. Его можно составить, потребовав, чтобы на приближенном решении задачи реализовывался минимум среднеквадратичной погрешность решения:

• 2

Ф(ХП, Х20, Х22) = ^ Х1{г, Х10, Х11, Х20> Х21), Х10, Х11, Х20 , Х21))] йг .

(20)

4-1

Таким образом, под приближенным решением уравнения (17) на отрезке а < г < Ь с начальным условием Х(г0 ) = Х0 понимается суперпозиция функций .X]. (г, Х10, Х11, Х20, Х21),

1 = 1, 2, ..., построенных в результате следующей процедуры: 1. На интервале а < г < Ь строим сетку а = г0 < г1 < г2 < ... < гн = Ь . Количество узлов сетки может меняться в ходе работы алгоритма.

г=г

г=г

1-1

Серия 1. Наземные транспортные средства, энергетические установки и двигатели.

2. Полагаем х10 = х0, / = 1.

3. Вычисляем длину интервала к. = — I ._х и строим функцию X. х10, х11, х20, х21), определенную формулой (17).

4. Вычисляем х11, х20 И х21, решая следующую задачу нелинейного программирования: минимизировать среднеквадратичную погрешность (20) при условиях

У х/ х10, х11, х20 , х21 )' х/ х10 , х11, х20 ' х21 )) _ 0 I

_ ^ \. (21) У (' х/' (' х10-х11' х20 , х21 )' х/ , х10 , х11, х20 , х21 )) _ 0 ]

5. Вычисляем меру относительной погрешности

о _ 2

8/ = Ъ/ ^ х/ и, х10 I х11, х20 ' х21 )' х/ х10, х11, х20 1 х21 ))] ' ^ (22)

•j'-l

И сравниваем С И Ômax . Если < Sj < Ômax , ТО полагаем j = j + 1, = Х20 и пеРеходим к п. 3, не меняя шага. Если ôj > Smax, то возвращаемся к п. 4 с уменьшенной длиной шага hj . Если же ô j < ômin, то полагаем / = J +1, x10 = x20 и переходим к п. 3, увеличив шаг интегрирования.

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

Обобщение изложенной методики на системы ОДУ сводится к замене скаляров x10, x11, x20 и x21 N-мерными векторами X10, Xn, X20 и X21, введением в подынтегральное

выражение функции (20) суммы квадратов левых частей рассматриваемой системы ОДУ, а также соответствующему увеличению числа ограничений в форме равенств.

Исследование изложенного алгоритма на стандартных тестовых задачах показало, что он обладает свойством абсолютной устойчивости: возмущения в начальных данных стремятся к нулю со скоростью (3/4) , где n - число шагов интегрирования. (Для неявного метода Эйлера, принимаемого при таком исследовании за эталон, этот показатель равен (12) ). К достоинствам разработанного алгоритма относится возможность контролировать

погрешность аппроксимации на каждом шаге интегрирования (см. п. 5), а также его способность отфильтровывать высокочастотные (нефизические) составляющие благодаря сглаживающему действию «целевой функции» (20).

Многозонные модели рабочих циклов бензиновых двигателей и дизелей содержат математические описания физических и химических процессов, зависящие от конструкции двигателя. Значительное число таких описаний (субмоделей) будут уточняться по мере развития теории ДВС. Существенно, что описанный алгоритм позволяет модернизировать и заменять эти субмодели (не затрагивая остальные - просто добавляя или исключая соответствующие ОДУ).

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

Литература

1. Poulos G., Heywood J.B. The Effect of Chamber Geometry on Spark-Ignition Engine Combustion. "SAE Techn. Pap. Ser.", 1983, 830334. pp. 24.

2. Ортега Дж., Пул У. Введение в численные методы решения дифференциальных уравне-

ний. М.: «Наука», 1986. 288 с. 3. Виноградов A.M., Красильщик И.С., Лычагин В.В. Введение в геометрию нелинейных дифференциальных уравнений. М.: «Наука», 1986. 336 с.

Применение координатных фотоприемников ФДК-124 в лазерных анализаторах загрязнений газов вредными примесями

Ащеульников Е.К., Климова Е.В., Тарасова Л.И.

Университет машиностроения 8(495) 223-05-23 (1587), а\г@татг.ти Аннотация. В статье рассматриваются вопросы применения лазерных источников света в анализаторах состава смесей выпускных газов автомобильных двигателей внутреннего сгорания, приведена функциональная схема установки, на которой возможно проводить эксперименты по испытанию и подбору квантовых генераторов для применения в анализаторах загрязнений газов примесями. Установлена возможность применения координатных фотоприемников типа ФДК-124 в анализаторах загрязнения газов микрочастицами.

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

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

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

Лучеиспускание лазера Ел во много раз больше лучеиспускания Еа обычного источника света, что видно из соотношения:

Е„ 1035 РО

ТТ =—Г°, (1)

Еа V

где: Р - выходная мощность излучения лазера;

О - добротность объемного оптического резонатора лазера на частоте V;

V - частота генерации лазера.

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

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

Например, сложно проводить отбор проб выбросов автомобильных двигателей автомо-