Статистические методы и идентификация математических моделей токсичности транспортных двигателей Текст научной статьи по специальности «Математика»

УДК 519.6:656.13

СТАТИСТИЧЕСКИЕ МЕТОДЫ И ИДЕНТИФИКАЦИЯ МАТЕМАТИЧЕСКИХ МОДЕЛЕЙ ТОКСИЧНОСТИ ТРАНСПОРТНЫХ ДВИГАТЕЛЕЙ

© 2005 г Е.Н. Ладоша, С.Н. Холодова, О.В. Яценко

Previous models of working process in piston engines are advanced in two ways: 1) intra-cylinder conditions are concerned with external ones, namely with the load, heat exchange and fuel supply as well as 2) an associative module to calculate the nitric oxides formation and destruction is developed. The model structuring is performed using asymptotic technique and its adjustment to known experimental data is realized statistically. As a result improved model enables really to predict numerically the detailed chemical composition of exhaust at arbitrary engine dynamics including non-stable operation. The accuracy of the model concerned with the estimation of ecological engine characteristics is much higher with respect to one of prototypes: the probable error respective to measured data is confined within 2-times region.

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

Кроме того, предыдущие попытки рассчитать токсичность поршневого ДВС при помощи компьютерных моделей [6, 7] лишь отчасти оправдали ожидания авторов: хорошо воспроизводя выбросы СО, углеводородов С„Нт, сажи и сернистых соединений, они давали на два порядка заниженные концентрации N0* в отработавших газах (ОГ). Поскольку в ранних версиях модели внутри-цилиндровых процессов «азотная» кинетика не учитывалась, неточность воспроизводства доли N0* в выхлопе связывалась с игнорированием химической неравновесности в параметрически однородном рабочем теле. Впоследствии [8] кинетическая схема была дополнена реакциями с участием азота, но это практически не повлияло на точность определения концентрации N0* в ОГ. В сложившейся ситуации, когда одно-трехзонные параметрически квазиоднородные модели термокинетики и химической термодинамики не позволяют адекватно воспроизводить токсичность ДВС по N0*, причи-

на их несовершенства может скрываться только в пространственном переосреднении [9] уравнений. Физически это означает, что замена истинного поля параметров в(г, /) средними по пространству величинами, например, <$(/)> или в(<г>, /) приводит к тому, что важные для целей моделирования связи уходят на подсеточный уровень, а само оно полностью утрачивает смысл. Избежать кардинального усложнения информационно-математической модели (ИММ) применительно ко всем компонентам триады «модель - алгоритм - программа» [10] удается в исключительных случаях: если в системе имеются (и выявлены) малые или большие параметры. При этом усовершенствование ИММ представляет собой добавление своего рода надстроек, сравнительно слабо связанных с базовой моделью, которая играет роль основания. В математике описанная техника эволюции моделей известна как метод сингулярных возмущений [11] - разновидность метода малого параметра [12].

В данной работе предлагается модификация модели [13] для учета 1) влияния переменной нагрузки и связанной с ней неравномерности вращения коленвала на процессы впрыска, распыливания, горения топлива и 2) особенностей кинетики образования окислов азота. Усовершенствования базируются на популярных в прикладной математике статистических методах. Ниже поочередно приводятся описания каждой из доработок модели [13].

Модель нестационарного режима работы ДВС

Основополагающие предположения усовершенствованной модели унаследованы от [13] и применительно к дизелям (с некоторыми оговорками также к инжекторным ДВС) формулируются следующим образом: 1) топливно-воздушный заряд в КС химически и термодинамически однороден; 2) впрыскиваемое топливо мо-нофракционно в соответствии с условной химической формулой С14Н30 (дизельное) или ОН^ (бензин); 3) реальный процесс горения адекватно описывается кинетической схемой, включающей возможные элементарные реакции между С-Н-0-содержащими веществами, имеющими в своем составе не более 3 атомов углерода, а также несколько брутто-реакций, аппроксимирующих термический распад (крекинг) тяжелых углеводородных фракций; 4) турбулентность корректно учитывается искусственным ограничением скорости элементарных химических процессов величиной темпа макро-микропереноса реагентов в зону и тепла - из зоны пламени; 5) испарение топливных капель происходит диффузионно; 6) теплообмен излучением незначителен. Сюда следует добавить известные предположения [2] касательно ки-

нематики кривошипно-шатунного механизма и оговорить характер сил, препятствующих вращению ДВС. После усовершенствований модель [13] принимает математическую форму:

ёП:

= /г (Пк , р, Т) + ^ (Г),

М

dT f ("к , Р, T){ / Cy - T}

dt

X "г

+ q(t)

dQ dA dQ

= q(t) _----^

dt dt dt dt

dQph

dA dt

sin Q t + f/2sin2 Q t

1 2 -+ 0,5-(1-cosQt+Y/2sin2 Qt

s-1

2 -

dQ

-QnRT, (1)

dt

■_aSw(T -Tw)

^L _ L F (t), ±П ^ dt f f dt dt

dmf 3Df

_ G(t)--J- mf,

dt

Л2

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

Усовершенствованная ИММ позволяет воспроизводить кинетику ВВ не только на установившихся режимах (результат такого расчета для локомотивного дизеля 10Д100 показан на рис. 1), но также при изменении нагрузки ДВС, величины и закона топливоподачи.

Динамическая и фазовые кривые Р(у) для одноцилиндровой версии этого двигателя при увеличении мощности от холостого хода до номинальной показаны на рис. 2. Разгон ДВС в вычислительном эксперименте осуществлялся увеличением цикловой подачи топлива в 4 раза; закон Оу оставался прежним.

3Бгтг ) = 8,/-/Х ,

I ёУё2 + (к! + к2 \dydt) ёф/Ж +

+ Р(у) Бр г(у) = 0, (2)

Р(У = кв КТ(УКУ , (3)

где пi = с/р, молекул/г - концентрации компонентов; /(пк, р, Т) - скорость наработки или распада химических веществ в реакциях, включая лимитирующий процесс нехимической - транспортно-диффузионной природы [6] 02 ^ 02 акт, скорость которого зависит от характера течения газа в заряде и задается как эмпирическая функция Яе для класса двигателей; Яе - число Рей-нольдса; - внешние источники газообразного вещества (при впрыске); q - доля теплового потока, не связанная с химическими реакциями; т/ и Ц - соответственно масса жидкой фазы и теплота испарения; ц - молекулярный вес _/'-го газового компонента; I - соответственно мгновенное значение площади поверхности теплообмена и характерного линейного размера КС; у и е - отношение радиуса кривошипа к длине шатуна и геометрическая степень сжатия; а, ТК - усредненные значения коэффициента теплоотдачи, температура поверхности КС; Б/, Л = Л(О) и 0(() - коэффициент диффузии паров топлива в вакуум, средний диаметр капель при распыливании форсункой и массовая скорость впрыска топлива; 8/ - символ Кронекера; Я - универсальная газовая постоянная; I - приведенный момент инерции нагруженного ДВС; Р(у), Б и г(у) - соответственно мгновенное давление в КС, площадь поршня и имеющая размерность длины функция передачи поступательного усилия - во вращение (это слагаемое в уравнении (3) соответствует одноцилиндровому ДВС); у -угол поворота и О = ёу/ё/ - частота вращения коленва-ла; к\ и к2 - коэффициенты сопротивления вращению; кв- постоянная Больцмана; прочие обозначения - согласно принятым в теплотехнике. В кинетической схеме уравнения учтены около пятисот газовых реакций между несколькими десятками веществ [14].

1,0 -

0,5

[СО],%

;

11111111111111 м 11111111111111111111

-T8D -эо о ао 160

Угол поворота коленеала, град до ВМТ

1,0 -,

0.5 -

i i i i i i i i I i i

— 1"SÜ —90

Угол поворота коленеала, град до ВМТ

i i I i 11

90

Рис. 1. Временная динамика усредненных по объему КС параметров рабочего цикла дизеля 10Д100 на номинальном режиме: а - относительное количество СО; б - объемная доля углеводородов СхНу в КС; ВМТ - верхняя мертвая точка

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

Смена стационарного режима в условиях расчета, результаты которого отображены на рис. 2, осуществляется за 20-30 оборотов коленвала, что соответствует примерно 2 с ускоренного движения городского автобуса. Разгон ДВС выражается увеличением примерно

а

вдвое каждой из величин - ёу/А и <Р(у) Sv г(у)> -среднего за период значения силового слагаемого в (2). Хотя конечный номинальный режим работы двигателя экологически чище (если, как это делается в данной работе, нормировать выброс ВВ на единицу израсходованного топлива), чем стартовый холостой ход, разгон двигателя сопровождается немонотонной динамикой выброса ВВ (в нашем случае окиси углерода). Концентрация СО в ОГ пропорциональна высоте плато за пиками на рис. 2 а и размеру «хвостика», отвечающего выпуску ОГ, на рис. 2 б. Максимальных значений этот показатель достигает «посреди» переходного процесса.

Приведенные на рис. 2 результаты моделирования относятся к простейшей ситуации - рассматривается одноцилиндровый двухтактный дизель и учитывается одно ВВ.

г [СО] , % 1.4 -

1.2 1.0 0.8 0.6 0.4 0.2

0

\hhhi

\W\hhhhhhhr\!

ш

ж

W

Л AlV

о

2000

4000 6000 8000 Я>, град. ПКВ

dlCО] о/

Рис. 2. Динамика переходного процесса в одноцилиндровом двухтактном дизеле: а - кинетическая кривая [СО]( у); б -соответствующий фазовый портрет. Выпуск ОГ при у = 120 0 после ВМТ

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

ного по относительной токсичности суммарного выброса ВВ не обязательно сохранят монотонность.

Влияние первого фактора проиллюстрировано на рис. 3: рассмотрен одноцилиндровый четырехтактный дизель в процессе быстрого перехода с режима номинальной мощности на холостой ход. Для простоты выпуск ОГ в данной модели полагается мгновенным в нижней мертвой точке у = 180 град. после ВМТ - верхней мертвой точки.

120rP(t>' МПа

100

г d[Pit)SDrit)]/dt

номинальная мощность

/

F(i)Sp/-(f)

Рис. 3. Силовая динамика одноцилиндрового четырехтактного дизеля с турбонаддувом типа тепловозного 10Д100: а - давление в КС, б - работа сжатия-расширения, в - фазовый портрет мощности

Модель образования окислов азота

Множественные экспериментальные и теоретические данные свидетельствуют, что определяющий брут-то-токсичность двигателя по N0* оксид азота N0 образуется исключительно в пламенной зоне горящего рабочего тела. В каждый момент активного горения объем

а

в

этой зоны сравнительно невелик (доли процента от полного объема, занимаемого рабочим телом), зато температура в ней приблизительно на 2200 К превышает среднюю по КС. Вследствие чрезвычайно резкой температурной зависимости реакция образования N0 протекает с сильным преимуществом в пламенной зоне: несмотря на ничтожный относительный объем последней, именно в ней происходит наработка всего N0. Кинетика наработки N0 в пламенной зоне ДВС с погрешностью ± 30 % описывается предложенным Я.Б. Зельдовичем [15] уравнением

-38000

= A • Р(у)[0](у)[Ы2](у)е т(у)+т™ -О ^ , Л ёу

(4)

где Т(у) - средняя температура рабочего тела; Q- относительное тепловыделение; А - константа, размерность и значение которой зависят от выбора единиц измерения концентраций реагентов; 71ппам- температура пламени. Это уравнение можно использовать для грубых оценок токсичности поршневых ДВС по NО.

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

-32000

Wp,

P(p)

T (p)

([NO]| равн. при T = T (p) - [NO](p))

рел О ^-Правы. ]

где п - порядок реакции гибели N0. Практически имеет смысл ограничиться двумя реалистичными ситуациями - гибелью окиси азота в моно- и бимолекулятных реакциях. Одновременно следует учесть, что равновесные концентрации N0 в условиях ДВС оказываются на два порядка ниже фактических, что позволяет считать [N0]^. при т = ту) ~ 0. Таким образом, получаем два сценария внутрицилиндровой брутто-кинетики N0. Согласно первому

-38000

, dQ

= A • P(p)[O](p)[N 2 ](p)e T(p)+T"

dt

•Q-

- B •

P(p)

-1

Q

-32000 T (p)

[NO].

dp

(5)

второму -

d[NO] = С • P(p)[O](p)[N 2](p)e T(P)+TnnaM • Q

dt dp

- D •

P(p)

-I

Q

-32000 T (p)

[NO]2.

(6)

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

Идентификация брутто-модели токсичности ДВС по окислам азота

Опишем процедуру определения параметров модели А, В, С и Б на примере уравнения (5). Предположим, что концентрация N0 в ОГ на момент их выпуска есть некоторый функционал Ф, зависящий от Р(у), Т(у), Q(у), [0](у), [N1 Р(у), ограничивающих значений уЛ]П, у^х и параметров А и В, т.е. Р^0]у = у^ = Ф{А, В, Р(у), Т(у), Q(у),

[0](у), [N2] Р(у), уть Ртах} . (7)

Теперь допустим, что сами функции и функционал Ф достаточно «хорошие» для того, чтобы (7) допускал приближенное представление [N0]^) = Ф{А, В, х}, где вектор х объединяет начальные условия расчета, а также параметры уП]П и утах. Последнее означает, что из рассмотрения исключаются неустойчивые или хаотические режимы работы ДВС [16], которые не представляют интереса для целей экологической оптимизации.

На следующем этапе формируем достаточно представительный вектор х (номенклатуру важных начальных параметров) и отбираем наиболее надежные и непротиворечивые среди имеющихся экспериментальных данных. Каждой реализации или конкретному набору значений компонентов х присваиваем номер ,, и интегрируем кинетические уравнения из [6,7] совместно с уравнением (5) для полигона параметров (А! ... АЛ) х (В! ... Вк). Наилучшими значениями А и В признаем те, которые дадут минимальное совокупное отклонение от всей рассматриваемой совокупности экспериментальных данных {х} по норме

2 [lgФ(А, B, Xj) - [NO](x; )]2 ^ min ,

(8)

i=1

Базой надстроечных моделей (5),(6) служат получаемые путем интегрирования уравнений из [6,7] эволюционные кривые Р(у), Т(у), Q(у), [0](у) и [К2](у). Небольшое число и отсутствие среди них кинетических кривых, отвечающих горению углеводородов, служат обоснованием применяемой нами техники наращивания модели. Фигурирующие в моделях (5),(6) константы А, В, С и Б являются параметрами, значения которых следует выбрать из условия наилучшей аппроксимации имеющихся в распоряжении авторов экспериментальных данных по N0*-токсичности ДВС. Чем больше первичных факторов или показателей используется при

где I - общее число опытных данных, используемых для инициализации модели. Выбор нормы логарифмического вида (8) обусловлен тем, что концентрации N0 вдвое большая и вдвое меньшая истинной следует считать одинаково плохими.

Практическое отыскание оптимальных значений А и В как решения системы уравнений

£ Ф{А, В, х}] - [N0],- | ■ Э Ф{А, В, х,} / Э 1&А = 0,

,=1

£ Ф{А, В, х,}] - [N0], | . д 1в Ф{А, В, х,} / д 1в В=0, (9)

,=1

затруднительно, поскольку не только требует полной переработки программного продукта [8], но и не допускает фрагментации соответствующей вычислительной процедуры. Отмеченное обстоятельство носит принципиальный характер и связано с незнанием характера межфункциональных связей, обозначенных как Ф{ А, В, х,}. Априорные сведения о подходящем шаблоне для

Ф{А, B, x,} позволили бы кардинально облегчить решение задачи (8). Выберем для аппроксимации истинных значений Ф{А, B, x,} шаблон Ф(А, B, x,) в виде линейной функциональной зависимости между логарифмами независимых и зависимой величин

lg Ф(А, B, x,) = a lgA + lg B . (10)

Представление (10) обосновано тем, что наработка и гибель основной части NO разнесены в пространстве и/или времени, следовательно, результирующий процесс можно рассматривать как их суперпозицию. Фигурирующие в уравнении (10) коэффициенты a и в служат своего рода проекцией вектора начальных и ограничивающих параметров x на параметрическое подпространство меньшей размерности, т.е. ux ^ u{a, в}, где стрелка символизирует параметрическое проектирование. Следующей гипотезой служит утверждение об информативности выбранного способа проектирования: отсутствие достаточных для этого оснований не следует непосредственно из постановки задачи и не может остаться незамеченным в процессе численных расчетов. С целью подтверждения или опровержения этой гипотезы ниже приводится анализ соответствующих данных.

Несмотря на формальное усложнение алгоритма: в нем по сравнению с методами непосредственного решения уравнений (9) появляются наборы дополнительных искомых коэффициентов a и в (по числу i опорных экспериментальных вектор-параметров x,), предлагаемая техника решения имеет ряд неоспоримых преимуществ. Во-первых, процедура решения расщепляется на два последовательных этапа: на первом для имеющихся параметрических полигонов ux, находятся проекции u{ a, в г}, а на втором по формулам

lg A =

í a,- lg[NO] i •í ßß -í ß, lg[NO]¡

i=1

i=1

i=1

•í a, ßr

i=1

ía •í ß2 -ííaßr

i=1

i=1

i=1

lg В =

I

порядков симметрично по обе стороны от этой точки (по каждому из направлений). Шаг дискретизации по ^ А, ^ В достаточно задать равным 1. Ограничиваясь шириной параметрического окна в 8 х 8 порядков, получаем изображенный на рис. 4 полигон узловых значений.

A4 B-4 A4 B4

Ao Во

A¡ Bk

A-A B-4 A-4 B4

í a • í ßt lg[NO]i - í a lg[N0] • í aß i

i=1 i=1 i=1 i=1

í a •íßß-\íaßi

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

Для вычисления агрегированных факторов влияния а и в строится вспомогательный полигон параметров \%А, х ^ Вк: удобно, выбрав базовое сочетание значений

А0, ^ В0}, очертить прямоугольные границы на 4-5

Рис. 4. Схема дискретизации вспомогательного проективного поля lg Aj x lg Bk

По результатам вычислительного эксперимента можно корректировать центровку параметрического окна, его ширину, а также параметрическое разрешение.

Чтобы определить величины а и Pi для /-го набора экспериментальных данных x, по рассчитанным путем интегрирования дифференциальных уравнений эволюционным кривым из [6-8] вычисляются значения Ф^ = ФА, Bk, x,), затем - оптимальные, т.е. минимизующие

функционал S [ lg Фк - а lgAj - Pi lg Bk ]2 значения а и

j к

Р.. Дифференцируя выражение функционала по переменным а и Pi и разрешая результирующую систему

линейных уравнений, получаем явные выражения

2

а= [ S lg Фк lgAj JS (lg Bk) -

j ,k k

- S lg Ф*' lg Bk ■ S lg Aj lg Bk] /D,

j,k j,k

p=[ S lg Фд lg Bk ■ K S (lgAj)2 -

j,k j

-S lg Ф*' lgAj S lg Aj lg Bk] /D,

j ,k j ,k

D = JK S (lgAj)2 S (lg Bk)2 - ( S lg Aj lg Bk)2] ,

j k j ,k

где J и K - соответственно число узловых ячеек (см. рис. 4) по А и B.

Поскольку предложенная процедура сопряжена с двукратным осреднением (в пределах полигона подгоночных параметров для ед иничного экспериментально -го набора x и по ансамблю таких наборов), она должна обладать высокой вычислительной устойчивостью. Побочным эффектом переосреднения является неспособность модели воспроизводить кинетику NO в малых областях параметрического пространства x, где чувствительность решения к влияющим факторам оказывается слишком высокой. Этим, в частности, объясняются

2

2

i=1

1=1

1=1

весьма жесткие требования к отбору экспериментальных данных х, для параметрической идентификации (можно сказать, инициализации) модели. На самом деле, если в совокупности инициализирующих векторов х, можно выделить две равномощные подгруппы, отражающие противоречивые или даже противоположные тренды, модель выявит не истинную ситуацию, а результат противоборства соответствующих статистик. Поэтому отбор сведений для инициализации данной модели следует возложить на экспертов по вопросам энергетической и транспортной экологии, а описанную процедуру осреднения модифицировать, введя статистический вес как фактор надежности /-го вектора инициализирующих данных х,. Техника взвешенного осреднения отличается от описанной лишь тем, что на заключительном этапе параметры А и В определяются как обеспечивающие минимум не (8), а модифицированного функционала «брутто-расхождения»

2 [ ¡я Ф(А, В, х) - [N0], ]2 м>,^ гат .

/=1

В процессе идентификации модели (5) необходимо вычислять а, и в, с точностью по меньшей мере до трех десятичных знаков: в противном случае сокращаются возможности для факторного анализа модели и результатов ВЭ. Очень полезен также контроль промежуточных результатов, наделенных определенным смыслом: в частности, величины а,, отвечающие за наработку N0 в ,-ом сценарии процесса, должны быть (если не все, то уж как правило) положительными, а агрегированные показатели гибели в г - отрицательными.

Вычисляя для каждого х, наряду с показателями а, и в, среднеквадратичное отклонение фактических данных ¡Я Ф. от их модельного представления а, ¡яА + в, ¡Я В*,

А , = [1/JK■ Z [lg Фк - a lgAj - ß■ lg Bk ] ]

2,1/2

j,k

(11)

Рис. 5. Погрешность функционального представления экспериментальных данных по МО-токсичности среднеоборотного транспортного дизеля [3] посредством шаблона (10)

можно оценить адекватность модели (10) истинному положению вещей. Кроме общей оценки качества упомянутой линейной аппроксимации посредством (11), интересно исследовать локализацию и характер имеющихся отклонений. Точечным аналогом (11) служит поле Д.1 = ¡я Ф. - а \%А] - в ¡Я Вк, вариант которого показан на рис. 5. Представленные здесь данные свидетельствуют, что использованный автором сравнительно скромный полигон экспериментальных данных х, [1-3, 5] позволяет настроить модель (5) для определения содержания N0 в ОГ среднеоборотных дизелей с погрешностью в пределах множителя, равного 2 (рис. 6).

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

Рис. 6. Сопоставление параметрического тренда [N0]^ = р^Ре) для дизеля Д-21А1, предсказанного в рамках данной модели, с неиспользованными при ее инициализации результатами экспериментальных исследований. Ширина коридора вероятных значений оценивается равной 4: от 1/2 до 2 [N0] р = ^пах- Сплошная линия - эксперимент [4]; □ - расчет данной работы; коридор вероятных значений ограничен штриховыми линиями и затемнен

Обозначая наработку N0 в наиболее общей форме как Ф(А, В, х,) и решая задачу обоими методами, мы можем составить два функционала

Дш = [1/1 2 [ ¡я Ф(А, В, х) -а ¡я А - вг ¡Я В ]2 ]1/2 ,

(12)

ANo = [1// Z [ a lgA + ß lg B- [NO], ]2 ]

2 1/2

i=1

(13)

Первый служит мерой того, насколько справедливо проектирование их ^ и{а, в соответствии с шаблоном (10) для всего выбранного параметрического полигона их,. Таким образом, формула (12) является обобщением (равнозначно осреднением) индивидуальной шаблонной погрешности (11) применительно ко всей совокупности опорных данных и{х,, [NO],}. Величину ДШ следует осмысливать как меру пригодности (10) в качестве функционального шаблона для lg Ф(А, B, х,) и в равной степени как меру эмпиричности «универсальных» подгоночных параметров A и B. Вычисляемое, согласно (13), значение Дда отвечает вероятным погрешностям выдаваемого идентифицированной моде-

i=1

лью прогноза касательно содержания окислов азота в ОГ. При этом, конечно, считается, что прогноз осуществляется для сочетания факторов, сходных с использованными для целей инициализации описанного N0*-блока. В этой величине «спрятана» ошибка воспроизводства частных результатов [N0], за счет введения универсальных констант А и В.

Агрегируя показатели ДШ и ДN0, получаем выражение для общей методической погрешности

Де = [1// 2 [ Ф(А В, х) - [N0], ]2]1/2 . (14)

1=1

Известный факт, что Де < ДШ + Дш имеет вполне понятную причину: введение всякого промежуточного этапа вычислений, особенно связанного с осреднением, неизбежно сопровождается накоплением ошибок. В данном методе, выбрав вид шаблона, мы вносим ошибку как в представление об индивидуальной ,-й реализации процесса (оценивается формулой (13)), так и в определяемые, согласно (12), совокупные сведения о пучке реализаций и{х,, [N0],}. Порождающее проблему противоречие заключается в том, что универсальная форма шаблона вводится как для частных реализаций, так и для их ансамбля. Как следствие подгоночные параметры оказываются переосредненными - аппроксимация нивелирует различия между отдельными реализациями. Таким образом, внося в модель элементы статистики, мы неизбежно наделяем ее специфическими недостатками. Для выбранного в целях апробации модели полигона и{х,, [N0],} погрешность аппроксимации Де примерно поровну распределилась между компонентами (12) и (13), что эквивалентно разумному (в соответствии с иерархией пространственно-временных характеристик элементарных процессов в системе) расщеплению задачи на слабо связанные фрагменты.

Отметим, что размерность фигурирующих в формулах (11)-(14) погрешностей аппроксимации Д,, Д^1, ДШ и др. при выбранном логарифмическом масштабе для коррелирующих переменных, характеризуемая как «порядок величины», представляется весьма удобной в инженерных оценках экологических показателей техники.

Располагая достаточно надежным и пространным массивом экспериментальных данных и{х,, [N0],}, все описанные действия целесообразно проделать также над моделью (1), четвертое уравнение. Затем, сравнивая соответствующие показатели качества аппроксимации -Д,, Д,ь ДШ, ... обеих моделей, можно отдать предпочтение тому или иному механизму гибели N0 на стадии догорания - расширения.

Описанный выше алгоритм инициализации ИММ образования окислов азота в ДВС можно классифицировать как метод каскадной (двухэтапной) факторизации.

Таким образом, усовершенствованная математическая модель рабочего цикла ДВС позволяет заметно уточнить оценки загрязнения городского воздушного бассейна плотными автотранспортными потоками без существенных экспериментальных и стендовых измерений, неизбежных при использовании действующих методов оценки [1-5]. Упрощение процедуры квалификации состояния наземных транспортных средств, а следовательно, и городского воздуха при помощи данной ИММ позволяет рассматривать ее как перспективное средство комплексного экологического мониторинга.

Литература

1. Якубовский Ю. Автомобильный транспорт и защита окружающей среды. М., 1979.

2. Колчин А.И., Демидов В.П. Расчет автомобильных и тракторных двигателей. М., 2002.

3. Беднарский В.В. Экологическая безопасность при эксплуатации и ремонте автомобилей. Ростов н/Д, 2003.

4. Марков В.А., Баширов Р.М., Габитов И.И. Токсичность отработавших газов дизелей. М., 2002.

5. Кульчицкий А.Р. Токсичность автомобильных и тракторных двигателей. М., 2004.

6. Булыгин Ю.И. и др. Расчет энергоэкологических параметров ДВС «ENGINE» / Свидетельство об официальной рег. ПрЭВМ № 2002610605. М.: РОСПАТЕНТ, 2002.

7. Жигулин И.Н., Ладоша Е.Н., Яценко О.В. Тепломассообмен в энергетических и транспортных системах. Ростов н/Д, 2002.

8. Яценко О.В. Информационно-математические модели внутрицилиндровых процессов в поршневом двигателе: механизмы, уравнения, термодинамические и кинетические коэффициенты / Свидетельство об официальной рег. БД № 2004620137. М.: РОСПАТЕНТ, 2004.

9. Колесниченко А.В., Маров М.Я. Турбулентность многокомпонентных сред М., 1999.

10. Самарский А.А., Михайлов А.П. Математическое моделирование: Идеи. Методы. Примеры. М., 2001.

11. Коул Дж. Методы возмущений в прикладной математике. М., 1972.

12. Тихонов А.Н. Система дифференциальных уравнений, содержащих малые параметры при производных // Мат. сб. 1952. Т. 31 (73). № 3. С. 575

13. Яценко О.В., Загороднюк В.Т. Компьютерное моделирование задач прикладной физико-химической динамики. Ростов н/Д, 2001.

14. Яценко О.В. Радиационные, кинетические и термодинамические константы для компьютерного моделирования взаимодействий технических систем с неравновесными газами рабочей и окружающей сред / Свидетельство об официальной рег. БД № 2003620238. М.: РОСПАТЕНТ, 2003.

15. Зельдович Я.Б., Садовников П.Я., Франк-Каменецкий Д.А. Окисление азота при горении. М.;Л., 1947.

16. Малинецкий Г.Г. Хаос. Структуры. Вычислительный эксперимент: Введение в нелинейную динамику. М., 2000.

Донской государственный технический университет_J2 нояб^ря 2004 г