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

Вычислительные технологии

Том 15, № 5, 2010

О погрешности измерения температуры с помощью термопары в случае больших градиентов*

А. Д. Рычков, В. Д. Лисейкин Институт вычислительных технологий СО РАН, Новосибирск, Россия Новосибирский государственный университет, Россия e-mail: rych@ict.nsc.ru, lvd@ict.nsc.ru

А. В. КОФАНОВ Новосибирский государственный университет, Россия e-mail: avkof 87@gmail. com

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

Ключевые слова: математическое моделирование, численные методы, теплопередача в твердых телах.

Введение

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

* Работа выполнена при поддержке междисциплинарного интеграционного проекта фундаментальных исследований СО РАН № 26, проекта Президиума РАН № 1.4, интеграционного проекта СО РАН совместно со сторонними научными организациями № 94 и РФФИ (грант № 09-01-12023).

1. Математическая модель взаимодействия тепловой волны с термопарой

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

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

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

— термоэлектрический эффект термопары определяется температурой, осредненной по объему этого шарового сегмента Ят;

— теплофизические свойства ТТ и материала термопары различны, постоянны и не зависят от температуры.

Процесс теплообмена термопары с ТТ рассматривается в области 0{0 < г < Яр, х3(Ь) < х < Ь, Ь > 0} (см, рис, 1), в которой термопара находится на оси симметрии. Здесь Яр, Ь — радиус внешней границы и правая граница области, расположенные достаточно далеко от головки термопары; х3(Ь) — абсцисса поверхности горения, которая перемещается вправо с постоянной скоростью ир так, что х8(Ь) = х8(0) + ирЬ. Начало координат находится в геометрическом центре шарового сегмента.

Уравнение теплопроводности в области О в цилиндрической системе координат записывается в дивергентной форме

д_

т

д дТ 1 д дТ

^ ^ ^ дх ^ дх ^ г дг дг

(1)

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

Т(г,х8(Ь), Ь) = Т3,

дТ (г, х, Ь)

дх

0,

дТ (г, х, Ь)

х=Ь

дг

дТ (г, х, Ь)

т=Я,

р

дг

0.

г=0

Рис. 1. Схема расположения термонары

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

Т(г,х, 0) = То + (Т - То) ехр^-ир(х - Х2(0))Сррр/\^ ,

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

То

2. Численный метод решения

Для численного решения уравнения (1) применялся метод конечных объемов, позволяющий проводить расчеты на произвольной конечно-разностной сетке. Уравнение записывалось в интегральной форме для произвольного фиксированного объема V:

д_

Ж

QdV - ф Е^ = 0,

(2)

V

где Q = СрТ, Е — поток тепла через ориентированный по нормали элемент dS поверхности ограничивающей объем V. Пусть в области О задана произвольная разностная сетка, каждая ячейка которой топологически эквивалентна прямоугольнику. Обозначим объем такой ячейки через V] и среднее значение величины ^ на н-м слое по времени, отнесенное к центру этой ячейки, — через Qnj•. Тогда уравнение (2) аппроксимируется следующим разностным соотношением со вторым порядком точности по времени и по пространству:

4Q 1 - зQnj+1

2т

У* -

(Е ■ ^ - (Е ■ ^

(Е ■ - (Е ■

0,

(3)

т

потоками тепла через соответствующие площади граней объема У] умноженные на единичные нормали к ним. Способ их вычисления описан в [4]. Полученная разностная схема нелинейна и для ее решения можно применить итерационную схему, основанную на введении псевдовремени на каждом слое по времени:

Qi,j - Qi

Л]

3Qj1,2 - 4Qnn? + Q:

,п— 1

1,]

т1

2т

У]

(Е ■ S)n++l1/2s - (Е ■ S)n_+

1 —1/2

(Е ■ s)n+l1/2 - (Е ■ s)n+l1/2

0,

(4)

где т1 — шаг то псевдовременн, в — номер итерации по псевдовремени. Для реализации (4) предлагается использовать схему расщепления по пространственным переменным [5] (нижние индексы опущены):

5s+1/2 - 5s

Vj - Aiis+1/2

- 4Q™, + Q;-1

Tl

%J

(F ■ S)n+i1/,s - (F ■ 1/2

\ra+1 ,s

, 3

2t

V,3 -

(F ■ s);+Î/2 - (F ■ S);-^'

^s+1 - ^s+1/2 T1

Vitj - À2^s+1, Q

,;+1,s+1

Q

;+1,s

5'

s + 1

(5)

Здесь 6е — поправки к величине ^ Ль Л2 — разностные операторы, учитывающие только вторые производные по соответствующим направлениям. После сходимости итераций по псевдовремени 6е = 0 и имеет место точная аппроксимация полного исходного уравнения. Граничные условия для 6е задаются следующим образом: на левой границе 6е = 0, на остальных — "мягкие" граничные условия; положения верхней и правой границ области выбирались таким образом, чтобы они не оказывали существенного влияния на распределение температуры в области, занятой термопарой.

3. Метод построения сетки

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

Основное уравнение, из которого находится искомое координатное преобразование s(£) : S2 — D (S2 — вычислительная область, D — физическая область), имеет вид

J 2 д

где

1 2 Т 7 ( ds%

s = х, s = р, J = aet -—т

R«M - « 92у - 9 « 92y J- « 92у

2 1У\ - 922 QÇ1QÇ1 912 QÇ1QÇ2 + 9ll flf 2fl£2 '

S . .

Функция Ф^(£)) в (6) обеспечивает возможность адаптации сетки,

D

D

бивалась на пять блоков, как показано на рис, 2, в каждом из которых сетка строилась отдельно. Такое разбиение обусловлено необходимостью согласования сетки с границей термопары. Первый блок (xs(0) < x < -Rm) представлял собой неравномерную прямоугольную сетку, В областях 1 и 2 сетка строилась с помощью обращенных уравнений (6) относительно управляющей метрики, отвечающей за сгущение к головке термопары, а в областях 3 и 4 — относительно метрики, отвечающей за сгущение к проволочке, соединяющейся с головкой термопары, В области 5 использовался одномерный алгоритм построения сетки.

При интегрировании уравнения (2) шаг по времени выбирался из условия "улавливания" поверхности горения в узел сетки, т. е, так, чтобы па каждом последующем шаге по времени положение левой границы х8(£ + Д£) = х8(£) + ирД£ совпадало с ближайшим к пей справа узлом разностной сетки. Тогда перестройка сетки в процессе решения пе требуется. Поэтому в области 5 (см, рис, 2) сетка строилась так, чтобы сеточные линии по направлению У были строго вертикальными. Степень сгущения выбиралась в соответствии с результатами решения задачи (1), Фрагмент построенной разностной сетки показан на рис, 3,

Разностная сетка имела размер 350 х 250 узлов в плоскости (х, г), что обеспечивало точность расчетов около 0,1%,

На каждом этапе расщепления разностной схемы (2) при вычислении величин ¿5+1/2, матрицы перехода являются диагональными, что позволило разработать параллельный алгоритм, основанный па одновременном вычислении с помощью скалярных прогонок указанных величин сначала по одному направлению (по оси ОХ), а затем по другому. При решении па кластерной вычислительной системе наиболее эффективным оказалось использование системы параллельных вычислений ОрепМР, которая па мпогоядерпых вычислительных узлах позволяет работать с общей памятью в пределах одного вычислительного узла. Это дает возможность избежать большого объема меж-

Рие. 3. Фрагмент разностной сетки

Рис. 4. Коэффициент ускорения на мнох'о-ядерном вычислительном узле: 1 идеальное, 2 реальное ускорение

процессорных обменов, характерных для системы параллельных вычислений MPI из-за

5

ления.

На рис, 4 (см, с, 103) приведена зависимость коэффициента реального ускорения (кривая 2) от числа ядер вычислительного узла, полученная на кластерной системе DL580G2 (ИВЦ Новосибирского государственного университета). Отклонение этого коэффициента от идеального его значения (прямая 1) связано с "накладными расходами", характерными для ОрепМР при создании и обслуживании параллельных ветвей. Следует заметить, что эти "накладные расходы" существенно зависят от соответствия структуры исполняемой пограммы и базовой иерархии памяти вычислительной системы и могут изменяться в достаточно широких пределах,

4. Некоторые результаты расчетов

Расчеты проводились для различных радиусов головки термопары Rm при следующих значениях основных параметров: радиус проволочки rm = 1.0 ■ 10-3ем, радиус внешней границы Rp = 5 ■ 10-2 см, начальное положение поверхности горения xs = -5 ■ 10-2ем, Использовались следующие теплофизичеекие характеристики сред: Cm = 0.2 кал/(г^К), Cp = 0.3кал/(г ■ К), Àm = 0.16кал/(см ■ с ■ К), Xp = 7.2 ■ 10-4кал/(см ■ с ■ К), pm = 8 г/см3, pp = 1.6 г/см3. Здесь индекс m относится к материалу термопары, индекс p — к материалу ТТ. Температура поверхности горения Ts = 650 К, начальная температура T0 = 300 К. Расчеты проводились до того момента времени, в который поверхность горения касалась головки термопары.

В качестве точки измерения температуры термопарой принимался геометрический центр ее сферической головки. Температура T(0, 0, t) определялась либо непосредственно в этой точке, либо путем осреднения температуры по всему объему головки. Различия температур, определенных обоими способами, не превышало 1 град. Температура на внешней границе, где отсутствовало влияние термопары, выбиралась в точке T(Rp, 0,t). На рис. 5 приведены изменения относительной разности указанных температур 5 = (T(0, 0, t) - T(Rp, 0, t)) ■ 100 % во времени для двух значений скорости горения и двух размеров головки термопары. Эта величина и является относительной погрешностью измерения температуры термопарой.

Причина появления погрешности связана с тем, что коэффициент теплопроводности материала термопары существенно выше коэффициента теплопроводности ТТ и "сток" тепла в проволочку от головки термопары изменяет распределение температуры в головке. Анализ полученных результатов для различных размеров головки термопары привел к достаточно неожиданному результату — погрешность измерения растет при уменьшении радиуса головки термопары. Последнее объясняется тем, что с уменьшением радиуса головки соответственно уменьшается отношение площади поверхности головки к площади сечения проволочки, т. е. возрастает интенсивность отвода тепла в проволочку. Поэтому естественное стремление уменьшить размер спая термопары вместо уменьшения погрешности измерения может приводить к ее росту.

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

температуры в ТТ при низких давлениях в окружающей его среде, когда скорость горения мала, могут быть существенно искажены. При высоких скоростях горения (при высоких давлениях) разность температур Т(0, 0,£) — Т(Ир, 0,£) становится отрицательной, поскольку в этом случае время воздействия тепловой волны па термопару мало и ее головка пе успевает прогреться.

Значительный интерес представляет изучение влияния термопары па величину скорости горения топлива, когда поверхность горения приближается к головке термопары, Поскольку скорость горения пропорциональна разности потоков тепла из газовой фазы к поверхности горения и от псе в топливо, то можно исследовать поведение величины Q = д(0, х8({), {)/д(Кр, х8({), ¿), оде д — соответствующие потоки тепла в топливо (см, штриховые линии па рис, 5), Из данных рис, 5 следует, что при приближении поверхности горения к головке термопары значение величины Q возрастает многократно,

5, %

12

VI

/

( 1 1 2 1

/ / / I 1 1

/

;

---" — --- __ /

0.2

0.4

0 <>'-

28

21

14

0.6 /, с

-0.5

-1.5

-2.5

■ 2 1

1 \ 1

\| Г 1

! / **

О

25

20

15

10

0.002 0.004 0.006 о

5

Рис. 5. Изменение погрешности измерения температуры при скорости горения ир = 0.1 см/с (а) и ир = 1см/с (б): Ят = 40мкм (1), 80мкм (2)

Изо- Г, К

линии

16 630

15 552.5

14 475

13 435.418

12 397.5

11 371.508

10 346.383

9 341.411

В 337.866

7 333.73

6 330.783

5 328.821

4 327.3

3 324.8

2 322.763

1 320

0.004 х, ОМ

Изо- Г. К

линии

16 625.418

15 613.935

14 600.836

13 576.255

12 551.673

11 540.968

10 527.091

9 516.995

8 505.889

7 503.912

6 502.509

5 500.938

4 499.228

3 493.58

2 485.016

1 453.345

-0.004

0.004 Л'см

Рис. 6. Изотермы для Кт = 40мкм; ир = 0.1 см/с (а), ир = 1см/с (б)

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

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

Проведенные исследования позволяют сделать следующие выводы:

1 — численное моделирование показало наличие существенной погрешности в измерении температуры в TT с помощью термопары, величина которой зависит как от скорости горения, так и от размера спая термопары;

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

Список литературы

[1] Зенин A.A. Об ошибках термопарных измерений пламен // Инж.-физ. журнал. 1962. Т. 5, № 5. С. 67-74.

[2] Зенин A.A. О теплообмене термопар в волне горения твердого топлива // ПМТФ. 1963. № 5. С. 125-131.

[31 Asay B.W., Son S.F., Dickson P.M. et al. An investigation of the dynamic response of thermocouples in inert and reacted condensed phase energetic materials // Propellants, Explosives, Pyrotechnics. 2005. Vol. 30, No. 3. P. 199-208.

[4] Vinokur M. An analysis of finite-difference and finite-volume formulations of conservation laws //J. Comput. Phvs. 1989. Vol. 81. P. 1-52.

[5] Яненко H.H. Метод дробных шагов решения многомерных задач математической физики. Новосибирск: Наука, 1967. 195 с.

[6] Liseikin V.D. A Computational Differential Geometry Approach to Grid Generation. Berlin: Springer, 2007.

Поступила в редакцию 28 ноября 2009 г.