Например, Бобцов

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

УДК 535.8
МОДЕЛИРОВАНИЕ ОПРЕДЕЛЕНИЯ ТЕРМОДИНАМИЧЕСКИХ ПАРАМЕТРОВ ВЫСОКОТЕМПЕРАТУРНОГО ГАЗОВОГО ОБЪЕМА ПАССИВНЫМ ДИСТАНЦИОННЫМ МЕТОДОМ

© 2010 г. А. В. Войцеховский, доктор физ.-мат. наук; О. К. Войцеховская, доктор физ.-мат.наук; Д. Е. Каширский; И. С. Суслова
Томский государственный университет, г. Томск
Е-mail: vok@elefot.tsu.ru
Статья посвящена созданию инженерной методики определения термодинамических параметров высокотемпературного однокомпонентного газа исходя из интенсивности излучения, поступившего на фотоприемник. Методика базируется на предварительном точном определении интенсивности излучения для конкретного спектрального диапазона и фиксированных интервалов температур и парциальных давлений. Последнее осуществлено теоретическим расчетом с помощью информационновычислительной системы на примере угарного газа.
Ключевые слова: газ, температура, парциальное давление, излучение.

Коды OCIS: 280.4991, 300.6390.

Поступила в редакцию 26.03.2010.

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

личных концентраций (10–6 – 1 атм) и температур (300 – 2500 K). Вторая – развитие методов решения прямой и обратной задач газоанализа для неоднородных многокомпонентных газовых сред в условиях значительной вариабельности концентраций смеси и широкого интервала температур с учетом характеристик высокочувствительных мультиспектральных фотоприемных устройств (ФПУ) для регистрации излучения в различных спектральных диапазонах. Такие устройства в настоящее время могут быть созданы, например, на основе полупроводниковых квантово-размерных структур.
Наибольший интерес в настоящее время представляет решение обратной задачи – газоанализа струи нагретых газов двигателя летательного аппарата по излучению, зарегистрированному в различных спектральных диапазонах мультиспектральным фотоприемником.
Эксперимент и теория показывают, что существует выраженная спектральная зависимость излучательных характеристик выхлопных газов, позволяющая осуществлять селективный прием и обработку характеристик излучения струи нагретых газов двигателя с целью определения состояния его функционирования, причем эта зависимость специфична для различных пространственных распределений температуры и парциальных давлений газовых компонентов

“Оптический журнал”, 77, 9, 2010

37

в выхлопе [1]. Проблема заключается в необходимости комплексного учета характеристик триады “источник–среда–приемник” при выборе наиболее эффективного алгоритма решения обратной задачи газоанализа. Оптимальный для рассматриваемой проблемы тип детекторов должен отличаться для каждого газового компонента смеси узкими спектральными областями фоточувствительности, которые могут смещаться в ИК диапазоне путем изменения свойств материала детектора. Следует принимать во внимание и быстродействие приемника. Наиболее простой вариант – неохлаждаемый тепловой приемник с набором оптических фильтров, обладающий инерционностью 10–2 – 10–3 с, в то время как ФПУ на квантовых структурах обеспечивает быстродействие 10–9 – 10–12 с [3].
Кроме того, постановка задачи накладывает жесткие граничные условия. С одной стороны, нижней границей по интенсивности приходящего на ФПУ излучения является пороговая мощность излучения ФПУ; с другой стороны, излучение абсолютно черного тела является верхней границей, при достижении которой информация об условиях среды теряется. Кроме того, исследования характеристик излучения нагретого газового облака и анализ переноса этого излучения в атмосфере должны выполняться с одновременным учетом эффектов молекулярного поглощения и рассеяния на аэрозольных сферических частицах с оценкой точности регистрации излучения отдельных компонентов смеси мультиспектральным ФПУ.
Результаты теоретических и экспериментальных исследований подтверждают адекватность представления полупроводниковых структур с квантовыми точками как полупроводника с промежуточной зоной, что обеспечивает одновременную регистрацию приходящего излучения в различных спектральных интервалах. В то же время применение наноструктур, содержащих квантовые точки в квантовых ямах, позволяет реализовать спектральные характеристики ФПУ в виде узких спектральных полос, что обеспечивает необходимую селективность приема [3, 4]. Таким образом, только выполнение мультидисциплинарных работ открывает возможность дистанционной диагностики функционирования двигателя, позволяющей судить о его состоянии и температурном режиме на основе анализа газового состава факела двигателя по спектральной зависимости его излучения для обеспечения контроля нормального режима функционирования двигателя.

Следовательно, методология дистанционной диагностики функционирования двигательных аппаратов по анализу их спектров излучения для предотвращения техногенных катастроф является недостаточно разработанной. Приближенные расчеты спектральных характеристик высокотемпературных газов развиваются в течение длительного времени и в основном представляют собой различные аппроксимации, связанные с упрощением спектральных зависимостей коэффициентов поглощения, функций пропускания или излучения основных газовых продуктов сгорания топлива (СО, СО2, H2O, NO) [5–7].
В работе определяется зависимость интенсивности излучения (ИИ) от термодинамических параметров источника излучения – высокотемпературного газового объема. Расчет ИИ необходимо проводить с учетом аппаратной функции ФПУ. Решение этой задачи невозможно без информационно-вычислительной поддержки – программно-вычислительных комплексов, ориентированных на дистанционный анализ функционирования двигателей по эмиссионным спектрам их выхлопов или факелов [8–11].
Точный расчет интенсивности излучения нагретого газового облака
В газовых средах с существенными градиентами парциальных давлений газов и температуры имеют место процессы молекулярного поглощения, рассеяния и переизлучения. Ситуация усложняется необходимостью учета прохождения излучения в атмосфере, что приводит к решению уравнения переноса, т. е. неоднородного линейного дифференциального уравнения с переменными коэффициентами.
Решение этого уравнения в общем виде известно, однако обратная задача – определение характеристик ядер этого уравнения – требует разработки соответствующего математического аппарата. Современные идеологии используют аналитические и численные подходы для построения решений, методы решения некорректных задач, а также различные асимптотики [12–14].
В данной работе точный расчет выполняется авторами на основе разрабатываемых уникальных информационных систем, генерирующих базы данных по параметрам спектральных линий основных газовых продуктов сгорания топлива для широкого интервала температур и давлений и осуществляющих прямой расчет (line by line) спектральных излучательных и

38 “Оптический журнал”, 77, 9, 2010

поглощательных характеристик газов (СО, СО2, H2O, NO) [8–11].
Разнообразие в терминологии приводит к

некоторым неясностям в определениях абсолют-

ных величин излучательных характеристик

реальных сред. В данной работе используются

основные определения в формулировках моно-

графии [15]. При высоких температурах (по-

рядка 1000 K) газовый объем обладает сущест-

венным собственным излучением, превышаю-

щим фоновое атмосферное. Интенсивность этого

излучения, согласно [15], может быть найдена

следующим образом:

ν2
I(ν,T) = ∫ B(ν,T)εν (T)dν,

(1)

ν1
где B(ν,T) – функция Планка, εν (T) – излуча-
тельная способность нагретого газового облака.

Соотношение (1) определяет интегральное зна-

чение приходящего на приемник излучения в

спектральном диапазоне ν1–ν2. Функция Планка для однородного газового

объема с температурой Т имеет вид

B(ν,T) = 2πhν3c2/(exp(hν/kT) −1), (2)

где h – постоянная Планка, ν – частота излучения, k – постоянная Больцмана. Размерность В(ν, Т) – [эрг/(см2 с см–1)], что соответствует определению спектральной энергетической светимости в [15] и излучательной плотности потока в [16] и используется в данной работе. В дальнейшем величину, рассчитанную по формуле (2), определим для простоты как “интенсивность излучения абсолютно черного тела” (ИИ АЧТ), а величину, определяемую формулой (1), – как “интенсивность излучения” (ИИ). Переход к размерности [(Вт/см2)] см–1] требует деления на 107 и деления на 2π при переходе к излучению в единичном телесном угле. Отметим, что в данной работе используются энергетические величины.
Понятие излучательная способность в приближении локального термодинамического равновесия определяется формулой

ε(ν,T) =1− exp[−k(ν,T)ρL],

(3)

где k(ν, T) – коэффициент поглощения газа [см–1атм–1], Т – температура [K], ρ – парциальное давление исследуемого газа [атм], L – длина трассы [см].
Уточним принятые упрощения. Формула (2) справедлива для случая однородной среды, когда парциальное давление газа и температура постоянны на всем протяжении оптической

трассы, и именно такой случай рассматривается в данной работе. В качестве параметров среды задаются температура, парциальное давление газа и длина оптической трассы. Излучающий слой состоит из молекул оксида углерода (СО), излучение направлено по нормали к приемному устройству. Угарный газ выбран как обладающий наиболее простым и хорошо изученным спектром. Поскольку рассмотрение ведется по замкнутой схеме (рассчитанные с заданными значениями температуры и парциальных давлений в прямой задаче характеристики поглощения и излучения СО используются как входные данные в обратной задаче), аргументированный анализ применимости методики обеспечен.
Выбор рабочих спектральных диапазонов
Для каждого конкретного газа необходимо выделить спектральный диапазон с максимальным значением прошедшего через атмосферу Земли излучения. Основные газовые компоненты факела двигателя (СО, СО2, H2O, NO) различаются по строению и относятся к линейным и асимметричным молекулам. Наиболее просты в расчетах спектры СО, поэтому моделирование было проведено для окиси углерода. Эта молекула обладает интенсивными излучательными полосами в областях 4,8 и 2,3 мкм.
Однако выбор спектрального диапазона не ограничивается только критерием величины прошедшего излучения. Необходимо учитывать характеристики приемного устройства, осуществляющего селективное детектирование сигнала. Перспективное ФПУ может быть реализовано в двух вариантах, различающихся конструкционными и материальными подходами. Первый вариант – в виде многослойной структуры со слоями различного состава, обеспечивающими фоточувствительность в требуемых спектральных интервалах, причем предыдущий слой является оптическим окном для последующего слоя. Второй – в виде квантово-размерных наноструктур, реализующих селективность приема за счет межподзонных оптических переходов.
Расчеты проводились для второго варианта ФПУ, который в свою очередь накладывает ограничения на спектральный диапазон регистрируемого излучения в пределах спектральной разрешающей способности фотоприемника. Если спектральный интервал характеризуется выражением λ0 ± Δλ, где λ0 – длина волны максимума ИИ, то это ограничение составляет

“Оптический журнал”, 77, 9, 2010

39

I, эрг см–2 с–1

3
150
2
100
1
50

0 500

1000

1500

2000 T, K

“Гладкая” зависимость ИИ от температуры для нескольких парциальных давлений. 1 – 0,5, 2 – 0,75, 3 – 1,0 атм.

Δλ/λ0 = 9 – 11% [3]. Поэтому для проведения расчетов выбраны два спектральных диапазона: (2180 ± 120) см–1 и (4300 ± 240) см–1. Для этих диапазонов в широком интервале температур Т (300–2100 K с шагом 100 K) и парциальных давлений газа ρ (10–4–1 атм) рассчитаны более 300 значений I(Δν, Т) при разных ρ и Т. Этот массив чисел использовался в качестве исходных данных для получения аппроксимационных соотношений зависимости интенсивности излучения от термодинамических параметров газа.
На рисунке представлена “гладкая” зависимость ИИ от температуры для нескольких парциальных давлений, что позволяет искать ее аппроксимацию полиномом.
Установление аппроксимационных соотношений для интенсивности излучения газов в зависимости
от термодинамических параметров среды
Аппроксимация полученных значений интенсивности излучения от термодинамических параметров газа может быть ориентирована на две задачи. Для первой достаточно найти такую формулу, которая бы наиболее точно (с погрешностью до 5%, поскольку средняя относительная квадратичная ошибка измерения ИИ составляет ±5% при доверительной вероятности 0,68 [17]) восстанавливала исходные значения ИИ во всем рассматриваемом интервале температур и парциальных давлений для быстрого расчета ИИ. Вторая требует учета ограничений, накладываемых необходимостью решения об-

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

∑ ∑I(ρ,T) = αiρi + βjTj ij

(4)

оказался не применим по следующим причинам. Различие в абсолютных значениях парциальных давлений ρ (10–4 – 1 атм) и температуры (300–2000 K) привело к развитию неустойчивости решения обратной задачи, поскольку диапазоны давлений и температур имеют слишком разновеликие значения, и определить константы модели, дающие приемлемые варианты восстановления (расхождение во всем интервале не более 5%) исходных значений ИИ I(ν), не удалось. Следующий вариант

∑I(ρ,T) = αi (ρT)i i

(5)

не принес успеха вследствие однозначности

I(ρ, T) для различных комбинаций парциальных

давлений и температуры в произведении ρТ.

В результате был предложен следующий

путь – метод последовательной аппроксимации

ИИ сначала по одному параметру (температу-

ре) при фиксации второго (парциальных давле-

ний), а затем аппроксимация полученных коэф-

фициентов функцией от второго параметра. При

этом наилучшее восстановление реализовалось

при аппроксимации интенсивности полиномами

пятой степени по обоим параметрам. Были по-

лучены 36 коэффициентов разложения, которые

и определяют зависимость ИИ одновременно

от двух параметров: парциальных давлений

(в интервале от 2,5×10–3 до 1 атм) и температуры

(в диапазоне 900–2100 K). Сужение интервалов

парциальных давлений и температуры вызвано

установлением пределов допустимой погреш-

ности восстановления 5%.

В итоге ИИ в указанных диапазонах опреде-

ляется формулой

∑ ∑I(ρ,

T)

=

5 i=0

⎝⎜⎛⎜⎜⎜j=50

αij

ρj

⎟⎟⎟⎟⎟⎠⎞Ti

(6)

с коэффициентами, значения которых приведены в табл. 1 и 2. Размерности коэффициентов αij соответствуют степеням ρ и Т.
По определенным коэффициентам полинома были восстановлены значения ИИ для заданных интервалов температур, качество восстановления исходных значений иллюстрируется табл. 3 и 4.

40 “Оптический журнал”, 77, 9, 2010

Таблица 1. Коэффициенты полинома αij для определения излучательной способности CO в спектральном интервале Δν = 2180 ± 120 см–1

i j
012345

0

–0,03246

–6,91578×10–4 3,79184×10–6 –4,08224×10–9 1,69601×10–12 –2,50723×10–16

1 1635,05011

–11,69589

0,02561

–1,71836×10–5 6,03406×10–9 –8,42903×10–13

2 15556,16824

–99,92358

0,20462

–1,5532×10–4 5,15116×10–8 –6,35584×10–12

3 –41203,90689 281,98757

–0,6313

5,4322×10–4 –2,05625×10–7 2,89912×10–11

4 41538,88295 –292,72802

0,67986

–6,11799×10–4 2,41316×10–7 –3,52721×10–11

5 –15164,50948 108,65257

–0,25734

2,36784×10–4 –9,52452×10–8 1,41536×10–11

Таблица 2. Коэффициенты полинома для определения излучательной способности CO в спектральном интервале Δν = 4300 ± 240 см–1

i j
012345

0 –5,00541×10–4

–5,35202

–6,10181

–5,91591

13,71317

–6,41881

1 3,08589×10–6

0,04201

0,04014

0,03055

–0,07921

0,03783

2 –6,35579×10–9 –1,1273×10–4 –8,58517×10–5 –4,77734×10–5 1,48064×10–4 –7,28544×10–5

3 5,51619×10–12 1,18931×10–7 6,53566×10–8 3,29347×10–4 –1,14979×10–7 5,81661×10–8

4 –2,13896×10–15 –4,06129×10–11 –9,05405×10–12 –1,91182×10–11 4,38292×10–11 –2,15535×10–11

5 3,06803×10–19 6,1252×10–15 –1,70913×10–15 4,2566×10–15 –6,49831×10–15 3,02307×10–15

Таблица 3. Анализ восстановления исходных значений ИИ для различных термодинамических параметров в полосе 4,8 мкм
Δν = 2180 ± 120 см–1

ρ, атм

Т = 900 K

Т = 1500 K

Т = 2100 K

εисх 0,0025 7,40562

εрасч 7,4241

δ, %

εисх

0,250146 20,0606

εрасч 20,1201

δ, %

εисх

0,296361 34,7558

εрасч

δ, %

34,7562 0,001219

0,01 29,8196 29,0754 2,495596 81,0468 80,4348 0,755076 139,793 139,3730 0,300451

0,05 158,383 157,7783 0,381788 427,16 426,4522 0,165698 718,98 717,9502 0,143228

0,25 962,515 959,4986 0,313387 2531,15 2533,1794 0,080176 3967,44 3966,6980 0,018703

1 3340,2 3346,4576 0,187341 10481,2 10505,3785 0,230685 16800,3 16824,3891 0,143385

Таблица 4. Анализ восстановления исходных значений ИИ для различных термодинамических параметров в полосе 2,3 мкм
Δν = 4300 ± 240 см–1

ρ, атм

Т = 900 K

Т = 1500 K

Т = 2100 K

εисх

εрасч

δ, %

0,0025 0,01177 0,0119778 1,7655

εисх 0,11583

εрасч 0,115895

δ, % 0,0564

εисх 0,36882

εрасч 0,36899

δ, % 0,0485

0,01 0,04727 0,0482138 1,9966 0,46502 0,465426 0,0874 1,47696 1,4805 0,0529

0,05 0,24032 0,2433621 1,2658 2,36886 2,37033

0,062 7,51173 7,515114 0,0451

0,25 1,23583 1,2411134 0,4275 12,7179 12,71914 0,0097 39,9159 39,92591 0,0251

1 4,99403 4,980402 0,2729 57,2835 57,26237 0,0369 178,288 178,457 0,0948

“Оптический журнал”, 77, 9, 2010

41

Максимальная ошибка в восстановлении этих значений составила для первого спектрального интервала 2,49% и минимальная – 0,0012%; для второго интервала – 2,03% и 0,005% соответственно.
Стандартное отклонение относительных ошибок определяется формулой
∑( )Δ = M (Iiâîññò − Iièñõ ) Iièñõ 2/(N − M) ⋅100%, (7) i=1
где N – размерность массива обрабатываемых значений ИИ (опорных точек), M – число определяемых констант полинома. Для спектральных интервалов (2180 ± 120) см–1 и (4300 ± 240) см–1 Δ равна 0,68% и 0,48% соответственно.
Погрешность предсказания значений ИИ для величин ρ и Т, попадающих в рассматриваемые интервалы, но не использованные в процедуре определения констант полинома (интерполяционных точек), не превышает 5%.

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

Полученные аппроксимационные соотно-

шения использованы для определения средних

значений парциальных давлений и температуры

из известных значений ИИ. Для каждой пары

значений ИИ (для двух диапазонов) решалась

система двух уравнений с целью определения

ρ и Т. Анализировалось восстановление этих

параметров в указанных диапазонах.

Так как в данной работе исследование про-

водится по замкнутой схеме и ИИ рассчитыва-

ются методом line-by-line, имеется возможность

проводить анализ достоверности решения об-

ратной задачи. Процедура определения ρ и Т

сводится к решению системы нелинейных алге-

браических уравнений (СНАУ)

⎪⎪⎨⎩⎪⎪⎧ff12((ρρ,,TT))

= =

I1, I2,

(8)

где f1(ρ,T) и f2(ρ,T) – функции вида (6), соответствующие вышеуказанным спектральным

интервалам.

В данной работе для поиска решения СНАУ

использовался наиболее распространенный

алгоритм оптимизации – алгоритм Левенберга–

Марквардта (Levenberg–Marquardt Algorithm,

LMA) [18], который является комбинацией про-

стейшего градиентного метода и метода Гаусса–

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

F(x) = (f1(x) − I1 )2 +(f2 (x) − I2 )2 =
= r12(x) + r22 (x) = r(x) 2,

(9)

где х = (ρ, T) – вектор искомых, r(x) = (r1(x), r2(x)) – вектор невязки. Вычисление х на очередном шаге осуществляется в соответствии с
выражением

xi+1= xi −(H + λI)−1∇r(xi ),

(10)

где Н – матрица Гессе, вычисленная в точке xi; λ – параметр, I – единичная матрица, ∇r(xi) – градиент функции.
Была решена обратная задача для опорных точек (значений ИИ, по которым строились аппроксимационные формулы). Анализ точности определения ρ и T приводится в табл. 5. Минимальная и максимальная относительные погрешности составили соответственно 0,0046 и 6,96 для ρ и 0,00047 и 1,43%для T. Среднеквадратическое отклонение абсолютной погрешности составило 0,0078 для ρ и 3,9907 для Т.
Также была решена обратная задача для интерполяционных точек. Анализ точности определения ρ и Т приводится в табл. 6. Минимальная и максимальная относительные погрешности для ρ составляют 0,0018 и 8,0345%, а для Т – 0,00041 и 1,6022% соответственно. Среднеквадратическое отклонение абсолютной погрешности для ρ составило 0,01498, а для T – 13,0279. Как видно из табл. 6, произошло увеличение относительных погрешностей. Так, если для парциальных давлений в первом наборе 60% точек соответствовали значения относительной

Таблица 5. Анализ точности определения ρ и Т для опорных точек (156 точек)

Интервал относительных погрешностей (δ, %)

ρ

Т

0 < δ < 0,01

2 (1,3%) 3 (1,9%)

0,01 < δ < 0,1

13 (8,3%) 55 (35,2%)

0,1 < δ < 1

97 (62,2%) 94 (60,3%)

1