CC BY
16
УДК 681.3:536.24.08
Решение обратной задачи теплопроводности
для исследования коротких линейных тепловых
Канд. физ.-мат. наук А. В. СЕРЯКОВ1, А. П. АЛЕКСЕЕВ
1 seryakovav@yandex.ru ООО «Рудетранссервис»
Представлены результаты исследований путем решения обратной задачи теплопроводности теплофизических характеристик (теплоемкости) испарителя С<п, коротких линейных тепловых труб (ТТ) с паровым каналом, подобным соплу Лаваля, и предназначенных для охлаждения космических аппаратов и спутников со строгим регулированием взлетной массы. Математическая постановка коэффициентной обратной задачи теплопроводности в одномерной системе координат дополняется результатами измерений температуры, находящейся в адиабатических условиях поверхности ТТ вдоль образующей во всем диапазоне температурных нагрузок, теплового сопротивления RTT, тепловой мощности испарителя Qev и передаваемой в вихревой проточный калориметр тепловой мощности конденсации Qcond при монотонном и близком к линейному во времени нагреве испарителя. При большом температурном напоре и начале кипения жидкого рабочего тела измеряемая температура внешней поверхности капиллярно-пористого испарителя близка к постоянной, и проводя решение коэффициентной обратной задачи теплопроводности с экспериментальными значениями температуры испарителя, удается получить численную оценку экстремальной теплоемкости работающего испарителя ТТ и удельную теплоту испарения кипящего рабочего тела и сравнить ее с табличными значениями.
Ключевые слова: короткие линейные тепловые трубы, обратная задача теплопроводности, монотонный метод нагрева, тепловое сопротивление, теплоемкость.
Информация о статье:
Поступила в редакцию 09.11.2021, принята к печати 17.01.2022 DOI: 10.17586/1606-4313-2022-21-1-83-97 Язык статьи — русский Для цитирования:
Серяков А. В., Алексеев А. П. Решение обратной задачи теплопроводности для исследования коротких линейных тепловых труб. // Вестник Международной академии холода. 2022. № 1. С. 83-97. DOI: 10.17586/1606-4313-2022-21-1-83-97
The article presents the results of the researching thermophysical characteristics (heat-absorbing capacity) for an evaporator Cv of short liner heat pipes (HPs) with vapor channel (in the form of the de Laval nozzle) by a solution of inverse heat conduction problem. The pipes in question are intended for cooling spacecrafts and satellites with strict takeoff mass regulation. Mathematical statement of the coefficient inverse heat conduction problem in one-dimensional coordinate system is supplemented by the results of temperature measurements for the surface of the HPs under adiabatic conditions along generator line across the entire range of the temperature loads, as well as the measurements of thermal resistance RTT heat capacity of the evaporator Qv and condensation heat capacity Qcond transferred to vortex continuous-flow calorimeter at monotonous and nearly liner in time heating of the evaporator. At great temperature drop and the beginning of working fluid boiling, the measured temperature of the external surface of the capillary porous evaporator is close to constant. Therefore, solving the coefficient inverse heat conduction problem with experimental values of the evaporator temperature we are able to obtain both an estimation for extreme heat-absorbing capacity of the HPs' evaporator under operation and the values for specific evaporation heat of the boiling working fluid and compare them with the tabulated values.
A solution of inverse heat conduction problem for investigating short liner heat pipes
Ph. D. A. V. SERYAKOV1, A. P. ALEKSEEV
1 seryakovav@yandex.ru
LLC Rudetransservice
Keywords: short linear heat pipes, inverse heat conduction problem, monotonous heating method, thermal resistance, heat-absorbing capacity.
Article info:
Received 09/11/2021, accepted 17/01/2022 DOI: 10.17586/1606-4313-2022-21-1-83-97 Article in Russian For citation:
Seryakov A. V, Alekseev A. P. A solution of inverse heat conduction problem for investigating short liner heat pipes. Journal of International Academy of Refrigeration. 2022. No 1. p. 83-97. DOI: 10.17586/1606-4313-2022-21-1-83-97
Введение
Современный уровень развития вычислительной техники и вычислительных технологий позволяет нам значительно расширить класс решаемых прикладных задач. Таким образом, научно-технические проблемы теплопередачи ТТ, которые традиционно пытались рассматривать аналитически, все чаще анализируются и решаются с помощью численных методов с использованием специализированного программного обеспечения для инженерных расчетов. Речь идет о применении вычислительных программ для решения обратных задач теплопроводности твердых тел в одномерной постановке и использовании полученных методов для анализа работы коротких ТТ с паровым каналом, подобным соплу Лаваля, и большим количеством рабочей жидкости в капиллярно-пористой вставке и испарителе.
Численные методы получили наибольшее применение благодаря ряду специфических преимуществ, главным из которых является их относительно простая реализация на компьютерах [1]-[3].
Обратные задачи характеризуются нестабильностью их решений, что проявляется в возникновении больших численных изменений в решениях при небольших изменениях исходных данных. Задачи такого типа называются неправильно поставленными задачами или некорректными задачами [4]-[8]. Большой раздел некорректных задач состоит из обратных задач, возникающих в случаях, когда необходимые начальные и граничные условия (или данные) для постановки прямой корректной задачи, например, прямой задачи теплопроводности, отсутствуют, но имеется некоторая дополнительная информация о решении, позволяющая сформулировать обратную задачу. К таким задачам относятся обратные задачи теплопроводности (ОЗТ) [9]-[11].
Общую нелинейную задачу нестационарной теплопроводности ТТ можно рассматривать как совокупность трех нелинейных задач: задача, в которой нелинейность возникает из-за температурной зависимости коэффициентов основного уравнения; задача, в которой граничные условия нелинейны, например, из-за испарения или кипения рабочей жидкости в испарителе; и задача, в которой нелинейность возникает из-за температурной зависимости внутренних источников тепла (стоков).
Границы между этими задачами весьма условны, так как некоторые преобразования позволяют переключаться с одного типа нелинейности на другой [12]-[13], что иногда упрощает использование разработанных методов для решения нелинейных задач.
Основной трудностью при приближенном решении некорректных задач является выбор параметров регуляризации для достижения устойчивости решения. Для его определения наиболее широко используются следующие
подходы [12]-[15]: выбор параметра регуляризации по остатку в разностном функционале; использование вариационного метода с вычислением функции Лагран-жа (Лагранжиана); итерационный метод, в котором параметром регуляризации является число итераций, соответствующее ошибке входных данных [15]. Эта последняя хорошо известная методология решения обратных задач теплопередачи достаточно эффективна и широко используется в космической технике, военной науке, авиации, но ее реализация требует быстродействующих компьютеров и очень большого объема вычислительной работы [12]-[13].
При решении обратных задач для уравнений математической физики широко используются градиентные итерационные методы в вариационной постановке обратной задачи [13]. В этой статье мы рассмотрим простейший градиентный итерационный метод для приближенного решения ретроспективной (с обратным временем) обратной задачи по расчету теплоемкости испарителя коротких ТТ с паровым каналом, аналогичным соплу Лаваля, с известным значением теплового сопротивления RTT или коэффициента теплопроводности Х^. Кроме того, в нашей обратной задаче начальное условие итерационно уточняется, т. е. на каждой итерации решается обычная краевая задача для уравнения теплопроводности с уточняемыми начальными условиями, что повышает точность и надежность расчетов.
Кроме того, тепловая мощность в испарителе ТТ зависит только от времени т, а температурное поле T (z) во всем объеме (геометрическом пространстве) ТТ внутри адиабатического калориметра экспериментально определяется с использованием поверхностных термопар медь-константан. Эти измерения гарантируют, что обратная задача имеет однозначное решение, но это решение нестабильно; поэтому для решения задачи необходима дополнительная регуляризация и использование вариационного метода [14]-[15] с множителями Лагран-жа для определения оптимальных значений параметров в разностном функционале. Этот метод дает быстро сходящиеся последовательные приближения точного решения.
Нами разработан и реализован в системе Fortran для ПК простой алгоритм численного решения одномерной обратной задачи теплопроводности для расчета теплоемкости испарителя малой мощности с кипящей рабочей жидкостью. В этом случае важным и контролируемым фактором является температура внешней поверхности капиллярно-пористого испарителя ТТ, которая близка к постоянному значению. При дальнейшем нагреве уровень кипящей рабочей жидкости в испарителе ТТ медленно снижается, и температура внешней поверхности также медленно снижается. И это позволяет нам рассчи-
тать экстремальное поведение теплоемкости испарителя ТТ и оценить удельную теплоту кипения рабочей жидкости в нем.
Кроме того, что очень важно, наш простой алгоритм удобен для его реализации в измерительной схеме адиабатического вакуумного калориметра для исследования коротких ТТ. Возможности этого алгоритма иллюстрируются решением нескольких конкретных задач для расчета теплоемкости коротких твердых тел при монотонном нагреве при умеренных температурах, в т. ч. модельную задачу с испарением диэтилового эфира и, следовательно, с переменными начальными условиями.
Экспериментальная часть
Все экспериментальные исследования тепловых характеристик проводились с помощью разработанных ранее и выполненных из нержавеющей стали коротких линейных ТТ с паровым каналом в форме сопла, близкого к соплу Лаваля, подробное описание которых неоднократно приводилось в предыдущих публикациях [16]-[18].
В качестве рабочей жидкости используется диэти-ловый эфир С4Н10О, имеющий температуру кипения при атмосферном давлении Тв=308,65 К (35,5 0С), температуру замерзания Тр= 156,95 К (-116,2 0С), и критические параметры Тс=466,55 К (193,4 0С), Рс=3,61 МПа. Объем пор капиллярно-пористой вставки определяется в процессе изготовления, и в нашем случае он равен 16,6210-6 м3, масса диэтилового эфира с плотностью 713,5 кг/м3 (20 °С) в порах вставки равна 11,858 10-3 кг. Ограничение в массе заполняющего ТТ диэтилового
эфира сделано таким образом, чтобы испаритель не был залит эфиром при высокой температурной нагрузке и начале кипения пленки на нижней поверхности крышки ТТ (поверхности испарителя).
Коэффициент заполнения ТТ, определяемый по отношению объема диэтилового эфира к общему объему ТТ, равен 16,6210-6 м3/3,1410-5 м3=0,529. Длина наших ТТ составляет £тт=0,1 м, диаметр _0=2Т0-2 м, при заполнении диэтиловым эфиром средняя плотность ТТ равна ртт = 1,871-103 кг/м3 и, вследствие равномерного распределения эфира при заправке ТТ, плотность испарителя близка к плотности всей ТТ целиком, ртт ~ реу.
Средняя изохорная теплоемкость ТТ при отсутствии химического взаимодействия между компонентами (металлом и диэтиловым эфиром) равна Стт = 1,15 кДж/(кг - К). Теплоемкость ТТ как целого является слабо меняющейся функцией температуры из-за равномерного распределения рабочей жидкости в зонах конденсации, транспортировки и испарения и близка к суммарной теплоемкости диэтилового эфира и нержавеющей стали с массовыми долями этих компонентов при небольших температурных напорах (вне области кипения диэтилового эфира).
Измерительный вакуумный адиабатический калориметр, совмещенный с вихревым проточным калориметром для исследований ТТ показан на рис. 1 [18].
Схема включает в себя обозначения: 1 — вихревой проточный калориметр; 2 — фланец крепления ТТ; 3 -стеклянная крышка; 4 — крепление крышки; 5 — опорная ТТ; 6 — плоский резистивный нагреватель; 7 — сливной штуцер калориметра; 8 — входной штуцер — за-
Рис. 1. Схема проведения измерений коэффициента теплопередачи, пульсационных и тепловых характеристик коротких ТТ в вакуумном адиабатическом калориметре, совмещенном с вихревым проточным калориметром Fig. 1. Scheme for the calculations of heat transfer coefficient, pulsation and thermal characteristics of short HPs in a vacuum adiabatic
calorimeter combined with vortex continuous-flow calorimeter
вихритель калориметра; 9 — герметичный ввод измерительных проводов; 10 — емкостные датчики для измерения толщины пленки жидкого конденсата; 11 — измерительный и эталонный генераторы, расположенные на верхних торцевых крышках ТТ; 12 — внешний цифровой генератор; 13 — усилитель мощности; 14 — цифровой осциллограф; 15 — компьютер; 16 — управляемый коммутатор; 17 — цифровой вольтметр; 18 — сосуд постоянного напора воды; 19 — генератор воздушных пузырьков; 20 — расходомер; 21 — сосуд Дьюара; 22 — вакуумная камера; 23 — адиабатическая оболочка; 24 — охранные нагреватели; 25 — дифференциальные термопары для измерения поверхностной температуры ТТ; 26 — дифференциальные термопары адиабатической системы ТТ.
Внешние теплообменные поверхности зоны конденсации измерительной ТТ снабжены изолированными термопарами и установлены на глубине 1 диаметра ТТ в вихревом проточном калориметре со стабильным потоком воды. Для обеспечения точного измерения тепловой мощности и увеличения коэффициента теплоотдачи ТТ струйный поток поступающей воды в калориметре закручивают, значения скорости и завихренности потока регистрируются с помощью пузырьков воздуха. Число Рейнольдса Яеса1 в калориметре с температурой воды Тса1=(293±0,03) К равно Recaг=3,68T03, число Нуссельта №са1=77,3, коэффициент теплообмена аса/= 2,4103 Вт/(м2К).
Испарители ТТ, также оснащенные термопарами из медь-константановой проволоки диаметром 0,05 мм, нагреваются с помощью плоского резистивного нагревателя Н2, и их температура поддерживается на величину 57, выше температуры кипения диэтилового эфира 308,65 К при атмосферном давлении. Температура нагревателя стабилизирована, и величина перегрева испарителей ТТ устанавливается в диапазоне 57=0-20 К, при этом тепловая мощность измерительной ТТ не превышает 200 Вт.
Для уменьшения тепловых потерь при работе с ТТ их помещают в вакуумную камеру 22 из нержавеющей стали (10-3 торр), где они дополнительно окружены тонкостенным медным адиабатическим экраном 23, внутренняя поверхность которого покрыта слоем никеля, а на внешней поверхности размещены 4 секции защитных нагревателей 24.
Величина неадиабатичности вблизи середины ТТ
_2
не превышает ~ 2 10 _2К , в области резистивного нагревателя Н2 и испарителей ТТ неадиабатичность достигает ~ 10-1К.
Основная ТТ, называемая измерительной, заполнена диэтиловым эфиром, а эталонная, которая полностью идентична основной ТТ, заполнена осушенным воздухом под давлением 1 бар с температурой точки росы ниже 233,15 К (-40 °С). Коэффициент теплопередачи КТТ2 второй ТТ не превышает 0,15% от первой (измерительной) КТТ1 и не учитывается. Вторая ТТ, полностью идентичная первой, выполняет эталонную функцию при измерениях в дифференциальном режиме толщины пленки конденсата в первой ТТ [20]-[21].
Проведенные измерения толщины пленки конденсата на внутренней поверхности верхней крышки ТТ, с использованием емкостных датчиков и разработанных
высокочастотных генераторов [16]-[18], очень информативны и дают интересные результаты. Большая толщина пленки (и высокое тепловое сопротивление ТТ) при низкой температурной нагрузке и резкое уменьшение толщины пленки (и значительное снижение теплового сопротивления ТТ) с увеличением температурной нагрузки, по-видимому, связаны с изменением направления вращения тороидального парового вихря вблизи поверхности конденсации ТТ [20]-[21]. Все подробности опытного подтверждения изменения направления вращения парового вихря вблизи плоской поверхности конденсации внутри парового канала ТТ и зависимость толщины пленки конденсата от температурного напора на испаритель приведены в [16]-[18].
На рис. 2 показано температурное поле вертикально ориентированной ТТ при монотонном и линейном во времени нагреве с помощью расположенного снизу капиллярно-пористого испарителя: 1 — точки, соответствующие экспериментальным значения температуры поверхности ТТ Тшг с паровым каналом, выполненным в форме сопла Лаваля, К; 2 — сплошная кривая, расчетные значения температуры Т, К, в сопле Лаваля, сформированном паровом канале ТТ. В нижней части показана половина поперечного сечения парового канала вдоль продольной оси Oz.
На рис. 3 приведены экспериментальные значения температуры наружной поверхности испарителя ТТ с максимальным заполнением диэтиловым эфиром в начале кипения и при увеличенном масштабе. Отчетливо видно слабое падение температуры заполненного кипящим диэтиловым эфиром испарителя (< 1 К) и значительно более сильное (резкое) падение температуры паровой фазы эфира над испарителем внутри выполненного в виде сопла Лаваля парового канала, при этом экспериментальные значения поверхностной температуры ТТ в адиабатической оболочке калориметра близки к расчетным значениям температуры пара [16]-[18], подробности расчетов которого приведены в разделе «Распределение температуры внутри парового канала».
1. Уравнение теплового баланса испарителя
Наиболее общим и информативным нестационарным уравнением энергии для описания работы плоского испарителя ТТ в режиме монотонного нагрева является уравнение Фурье — Кирхгофа [5], которое содержит изохорную теплоемкость Сеу, Дж/ (кг-К), и среднюю плотность испарителя ре„, кг/м3 ; скорость потока паровой и жидкой фаз в векторной форме V, м/с, внутри капиллярно-пористого испарителя; диссипативные функции Ф^; изменение энергии при расширении движущегося пара Р^к (V); увеличение давления Ре„ в паровом канале; и без учета поправки на диффузию массы и импульса его записывают следующим образом:
Се^е^((+ V СТ^
= LeVF(())^ + div(V)ЧРеУ + Шу(1ттСТ) + е . (1)
Среднюю теплоемкость Сеу испарителя ТТ в случае химически не взаимодействующих компонентов оцени-
Рис. 2. Экспериментальные значения температуры поверхности ТТ вдоль образующей и расчетные значения температуры пара диэтилового эфира внутри образованного парового канала, подобного соплу Лаваля Fig. 2. Experimental values for the temperature of the HPs' surface along generator line and design values of diethyl ether vapor into a generated vapor channel in the de Laval nozzle form
335
a: 330
£ 325
320
, i i i
1*1 * ; 1
1 > 1 1
10
Длина, мм
15
20
Рис. 3. Экспериментальные значения температуры поверхности ТТ с максимально заполненным кипящим диэтиловым эфиром капиллярно-пористым испарителем в увеличенном масштабе Fig. 3. Scaled-up experimental values for the temperature of the HPs' surface with capillary porous evaporator maximum filled by boiling diethyl ether
вают обычным способом по формуле (2), где г безразмерная координата вдоль продольной оси ТТ:
Се, = Csc + хе,С,р + (1 - хе, ) С/еу + См>;
г . (2)
г =
L
HP
Qev = 1H 2 I
( T 2
и Э г
= EH 2 - CH 2 T ev =
z=0
(TH2 - Tfev )Fev (z)
= —\__uz-J ! ev\ ! = iig,4 gT .
d c + dwev
К к
fev 1 fev
Tfev = Th2-DTc-ATWeV- ATfev = 329,75 K(56,6 °C).
Qevdc
Qevd и
DTc =^ = 5,4 K ATWev = fFr = 18,8 K
^ЛЛ ¡-Lev
1 F
= Qevd fev
fev = F
fev ev
AT fev =
= (35...0,55 )K
(4)
При увеличении температурного напора в начальный период нагрева капиллярно-пористого испарителя ТТ смачивание сетчатого каркаса диэтиловым эфиром ухудшается, так как угол смачивания 9 увеличивается с повышением температуры. Этот эффект термочувствителен и уже при относительно небольшом нагреве & = Те, - Тв ~ (3 - 5) К смачивание ухудшается, уровень диэтилового эфира в сетчатом испарителе снижается и дальнейший нагрев ТТ & = Те, - Тв > 10 К через плоскую поверхность нижней крышки происходит в тонкий слой кипящего диэтилового эфира. Толщина слоя эфира уменьшается от величины ~ 3,5 мм в начале нагрева до (1 - 2) -10 -2 мм при максимально интенсивном кипении, при этом поток пара становится стационарным, и это значительно упрощает уравнение (1), и максимально достижимый результат в наших условиях выглядит следующим образом:
Уравнение тепловой энергии работающего испарителя ТТ при монотонном режиме нагрева в приближении малых теплопотерь кТТ (Тэт - Т!!Ие{) в адиабатическом калориметре может быть представлено с использованием энтальпии следующим образом: (=2 _ аеу (т) = Не, (т) + е т ФФ + кЭт (Тэт - ТзЫ). (5) (=1
Сама величина энтальпии Не, (т) при монотонном режиме нагрева испарителя также может быть представлена с использованием термодинамического уравнения:
НеV (т) = ОрТ(Тв) + Ор (т)С,р (Те, - Тс ) +
+о (т) Ср! (Те,-Те,). (6)
Выражение для тепловой энергии, выделяемой в области конденсации ТТ, записывается аналогично с помощью уравнения баланса тепла:
,dxe-
dTee
(3)
Qcond (t) = Gvpr (TB )ddv + GvpCvp (1 - Xev Y-jf +
+1 sc (1 - П)-2TгFev ()Lev - I VFjfl - 4T (( -TsM).
dz i=1
Проведенные ранее измерения и расчеты дают возможность оценить коэффициент теплопередачи в вихревом проточном калориметре асаг = 2,3103 Вт/(м2 • К),
&сопа - 5103 Вт/ (м2-К), и окончательное выражение для максимального значения теплоты конденсации в нашей ТТ:
Qcond
{(Г cond — Tc
cal ))cond
y fcond
hcal
cond
= 115,4 Вт
F {()
G = JL
vp 2 —X
M 2 pR
fcond
pvp =
. pvp
(8)
P_ P
RT' Pv
{k — 1)
1/k
= 1^1 ; P - P
fcond'
1 +
1
( 2 2 u — uv
k RTv
2 2 Ц v — v
k / {k —1)
(11)
Для решения сложного уравнения Фурье-Кирхгофа необходимо определить поле скоростей жидкой и паровой фаз в капиллярно-пористом испарителе с использованием уравнений Навье — Стокса, и провести стыковку с течением пара в паровом канале ТТ, это будет сделано позже.
Для упрощения всех последующих расчетов и вычислений среднюю температуру ТТ в поперечном сечении в вакуумном адиабатическом калориметре в любой точке на высоте ц = внутри парового канала мож-
но считать равной температуре поверхности ТТ, измеренной на этой высоте Т (,т) в вакуумной камере _ V Ьиг
адиабатического калориметра:
(ЛТтт = Ттт () - ТТТ ((). (9)
Паровой поток, обусловленный межфазным массо-переносом, может быть оценен в случае небольших отклонений от равновесия с помощью уравнения Герца — Кнудсена [22]-[23]:
ц1/2 Ж Р р^. Ц
(10)
где диэтиловый эфир С4Н100, выбранный в качестве рабочей жидкости, имеет молярную массу М=74,1216 г/моль; R — универсальная газовая постоянная; £ — коэффициент конденсации, определяемый как отношение числа конденсирующихся молекул пара, к общему числу молекул, которые попадают на поверхность пленки жидкости, при условии £ < 1.
Происходит перераспределение жидкой массы ди-этилового эфира внутри парового канала ТТ, при этом теплоемкость испарителя увеличивается за счет кипения медленно уменьшающегося по толщине слоя (в дальнейшем пленки) эфира. Близкая к постоянному значению температура Т^ слоя кипящего диэтилового эфира и зависимость от времени толщины слоя определяются с помощью внешних термопар, контролирующих температуру внешней поверхности ТТ внутри вакуумной камеры. Результаты измерений температуры с помощью термопар в момент начала кипения показаны на рис. 2 и рис. 3.
2. Распределение температуры внутри парового канала
Проведенные ранее расчеты скорости потока и плотности пара диэтилового эфира с использованием системы уравнений Навье-Стокса в коротком паровом канале ТТ, подобном соплу Лаваля [16, 18], позволяют нам рассчитать распределение температуры в этом паровом канале. Соотношение между давлением, плотностью и температурой конденсирующегося пара в первом приближении может быть задано уравнением состояния идеального газа в следующем виде [19]-[21]:
где k= 1,31 — значение адиабатического индекса пара диэтилового эфира.
Расчеты температуры движущегося пара в коротком паровом канале ТТ, выполненном в форме сопла Лаваля, были выполнены с помощью программного кода ANSYS\ CFD Fluent 6.3.26-20090623 [24] и с использованием уравнения (11) и других, приведенных в работах [19]-[21] более полных расчетных уравнений. С учетом табличных значений [25]-[27] давления пара диэтилового эфира p^," при коэффициенте влажности g dr = 0,2 было получено распределение температуры внутри парового канала вдоль продольной оси z при начале кипения в испарителе ТТ, как показано на рис. 2 и 3, при перегреве St = Tev — TB = 11 K:
f P #—1)k T=% Ц . (12)
Полученные расчетные значения температуры пара Tp вдоль продольной оси парового канала ТТ при температуре нижнего слоя диэтилового эфира в испарителе Tfev=329,75 К (56,6 °C) близки к экспериментальным значениям поверхностной температуры ТТ в адиабатической вакуумной камере. Погрешность расчета составляет 0,1%, погрешность измерения температуры поверхности ТТ составляет менее 0,3 K. Нагрев и расчеты проводились в диапазоне температур испарителя 298-348 К (25-75 °C).
Внешняя поверхность корпуса ТТ в области многослойного сетчатого испарителя высотой 3,5 мм, заполненного кипящим диэтиловым эфиром, характеризуется близкой к постоянной температуре при монотонном нагреве, что хорошо видно на рис. 2. При дальнейшем нагреве толщина слоя диэтилового эфира в испарителе и, соответственно, длина (высота) области постоянной температуры уменьшается.
Внутри парового канала температура резко снижается, падение температуры достигает ДТТ=35 К из-за охлаждения в области конденсации и влияния формы сопла Лаваля на течение пара в паровом канале ТТ. Результаты экспериментального анализа распределения температуры поверхности ТТ вдоль образующей также подтверждают нелинейный характер распределения температуры в зависимости от продольной координаты z.
3. Тепловое сопротивление ТТ
Экспериментальные значения в стационарном состоянии теплового сопротивления RTT короткой ТТ определяются как теплофизическая характеристика теплопе-редающего устройства в целом [28], при постоянной величине температурного напора на испарителе bt=Tev - TB и вычисляется с помощью формулы (13) и вакуумного адиабатического калориметра, в котором все параметры работающих ТТ, включая температуру и передаваемую тепловую мощность, измеряют с использованием термопар и калориметра, показанного на рис. 1.
Рис. 4. Тепловое сопротивление RTT, зависящее от перегрева испарителя St = Tev — TB относительно температуры кипения диэтилового эфира при атмосферном давлении: 1 — экспериментальные стационарные значения RTT коротких ТТ с паровым каналом, выполненным в форме сопла Лаваля; 2 — экспериментальные значения теплового сопротивления Rjy той же ТТ, полученные в режиме непрерывного нагрева с близкой к линейной во времени скоростью 3 10-3 K/c Fig. 4. Thermal resistance RTT depending on evaporator reheating St = Tev — TB in relation to diethyl ether boiling temperature at atmospheric pressure: 1 — experimental stationary values RTT of short HPs with vapor channel in the form of the de Laval nozzle; 2 — experimental values of thermal resistance RTT of the same HPs obtained with continuous heating at a rate of 310-3 K/s close to linear in time
RTT в режиме непрерывного нагрева, характеризуются направлением теплового потока, совпадающим с направлением теплового потока при измерениях во всех стационарных состояниях.
Полиномиальное уравнение (14), описывающее в безразмерной форме экспериментальные значения теплового сопротивления RTTk для момента времени тк короткой ТТ с паровым каналом, подобным соплу Лаваля, в зависимости от величины перегрева испарителя 5t: 0 < 5t < 20 (также в безразмерной форме), выглядит следующим образом:
nR
RTTk (dt) = ее Ri-ш (dtГ1 =-2,0795621 • 10 -8 (dt)7 + i=1
+1,6029662 • 10-6 (dt)6 -4,9921411 • 10-5 (dt)5 +
+8,0489929 • 10-4 (dt)4 - 0.0071936 • (dt)3 +
+0,03633406 • (dt)2 - 0.1113127 • dt + 0,2702057.
nR < 7. (14)
Стандартное отклонение с=0,0024929, критерий Фишера R2=0,9980243.
Полученные экспериментальные результаты теплового сопротивления при большом температурном напоре позволяют оценить характерную общую толщину пленки диэтилового эфира на поверхности испарителя (кипящей в испарителе) и поверхности конденсации коротких ТТ с паровым каналом, подобным соплу Лаваля, и плоскими верхними и нижними крышками следующим образом (без теплового сопротивления крышек):
" d cond ~ RTT1 fF =
= 4 10 -2—0,136 Вт
Вт
м2 • K
3,14• 10-4м2 J 1.7 • 10-6м'
Полная температурная нагрузка на испаритель 5t=Tev — Тв=0-20 К [16] (относительно температуры кипения эфира при атмосферном давлении) задается и поддерживается с помощью прецизионного регулятора температуры с шагом по температуре 5/=0,5 К, длительная нестабильность не превышает ±110-2 К.
Экспериментальные значения теплового сопротивления RTT, полученные при монотонном нагреве испарителя ТТ с близкой к линейной во времени скоростью 310-3 К/с с помощью специального прецизионного регулятора температуры и системы коммутаторов с компьютерным управлением, поверхностных термопар, измерения мощности нагревателя Н2 формируют массивы данных о температурах поверхности ТТ в вакуумном адиабатическом калориметре, поверхности конденсации, мощности нагревателя и тепловых потоках Qev (т) и Qcond (т) для решения уравнений теплопроводности и вычисления теплового сопротивления по формуле (13):
т _ т
2ет 1 СОПй
RTT (t) = -
Qei
Из анализа графика, показанного на рис. 4, следует, что толщина пленки диэтилового эфира на испарителе близка к постоянному значению в диапазоне тепловых нагрузок 5t =10-20 К, и при дальнейшем увеличении температурного напора 5t испаритель начинает высыхать и наблюдается резкое увеличение теплового сопротивления RTT.
Кроме того, когда испаритель ТТ монотонно нагревается и толщина слоя диэтилового эфира внутри испарителя 5fev (т) монотонно уменьшается, в диапазоне толщин 0 < z < Lev=0,035LTT, ограниченном продольным размером толщины испарителя Lev, функциональная зависимость теплового сопротивления RTT от координаты z: RTT=RTT (z) представляет практический интерес и имеет положительный знак:
-1
dRrr =3Rel ЗЭт Ц = Э%г ЖЭг Ч > 0.
Эг dt U dz J dt и at J '
(15)
(13)
Небольшой размер и низкие скорости нагрева гарантируют, что измерения RTT проводились в условиях температурной однородности и локального теплового равновесия. Полученные значения теплового сопротивления
0 J z J Lev = 0,035 Lll.
Распределение экспериментальных значений температуры поверхности ТТ, показанные на рис. 2 и 3, (dt/dz < 0) в начале кипения диэтилового эфира в насыщенном испарителе при монотонном нагреве, и уменьшающиеся значения теплового сопротивления dRTT/dt < 0 согласно экспериментальным данным, приведенным
на рис. 4, показывают положительное значение производной теплового сопротивления дRTT/дz > 0 по координате г в уравнении (15). И это подтверждает гипотезу о том, что зависимость теплового сопротивления RTT от продольной координаты г в интервале толщин испарителя г: 0 < г < Lev, характеризуется монотонным положительным ростом и единой функциональной зависимостью RTT=RTT (г). Это означает, что тепловое сопротивление RTT (г) может быть применено для решения уравнения теплопроводности во всем интервале толщин 0 < < Lev=3,5 мм капиллярно-пористого испарителя. То же самое можно сказать о теплоемкости испарителя.
Математическая модель обратной задачи теплопроводности
Математическая постановка обратной задачи восстановления (расчета) теплоемкости работающей ТТ может быть представлена следующим образом. Пусть имеется вертикально ориентированная ТТ длиной LTT= 100 мм и диаметром 20 мм, В^ / ^^=0,2, в вакуумном адиабатическом калориметре, температурное поле внутри ТТ рассматриваем в одномерной осесимметрич-ной постановке:
1 d 1
dt (z ,4)
z дzRev (() дZ
: Cev ((,tk);
z = -
-< 0,035.
LT
Для экспериментально определенных граничных условий обратная задача теплопроводности для уравнения (16) обычно некорректна в классическом смысле [9]-[11]. Чтобы привести ОЗТ к условно корректной формулировке, мы ограничиваем класс приемлемых решений набором кусочно-регулярных аппроксимирующих зависимостей и применяем пошаговый принцип естественной регуляризации для решения ОЗТ [6, 9].
Распределение температуры цилиндрического корпуса ТТ внутри вакуумной камеры с адиабатической оболочкой и вблизи плоского резистивного нагревателя симметрично:
dt(z,t) dt (z,t)
dx dy
z=0 z=0
= 0.
(18)
Начальные граничные условия в виде экспериментальных значений поверхностной температуры Ттт (zi)
/1 ^ у 'шг для решения уравнения (16) с учетом всех составляющих
теплового с°пр°тивления К^Лса1 ,
включая контактное сопротивление и ненулевую теплоемкость резистивного нагревателя СН2, (см. рис. 1), записаны стандартным способом [16, 18]:
* (0,Т* ) е = ТеV К) + (К + К + ) К ). (19)
(16)
Тепловое сопротивление испарителя Rev (испарительного фрагмента ТТ) является неотъемлемой частью теплового сопротивления RTT всей ТТ целиком, и величину теплового сопротивления испарителя при кипении (при температурном напоре 5T > 11 К) диэтилового эфира можно оценить, используя экспериментальные результаты, представленные на рис. 4. Учитывая близкие толщины пленок диэтилового эфира в испарителе и на поверхности конденсации ТТ при пленочном кипении §еу ~ dcond, тепловое сопротивление испарителя (испарительного фрагмента) составит при этом не более Rev ~ 0,02 К/Вт.
Потери энергии из-за трения в паровом канале и неабсолютная адиабатичность калориметра приводят к тому, что тепловая мощность Qev и Qcond отличаются друг от друга:
( - Qcond)) * 2,6% .
И тепловая мощность Qev становится функцией безразмерной вертикальной координаты z и времени при монотонном нагреве:
йе, = йе, ((,*). (17)
Для устойчивости решения обратной задачи теплопроводности (16) наиболее приемлемой ее формулировкой будет задание переопределенных граничных условий, с помощью которых и будет определено единство трех функций — температуры 1ШГ(,1к), теплоемкости Се, (,1к) и теплового сопротивления Ке, (,1к) (или теплопроводности 1ет (,1к), причем все функции будут аналитическими, например, представленными в виде полиномов или кусочно-гладких сплайн-функций.
t (0,0 ) = to;
t (0,tk ) = '0 +t k • 3 •10 -3К/с. ' (LTT ,xk) ='cond = Tcond (tk) + + (Rcal + Rw + Rfcond + Rcond ) Qcond (tk ).
(20)
(21)
Мы обозначаем величину теплового потока внутри ТТ в следующей стандартной форме:
9ev (Z,tk ) = Z^TT
(dt ((
dz
= EH2 - CH2Tev (z,Tk ) =
^z=0
(22)
= (2 - Tfev % (z) + V
0 fev 1
Тепловая мощность Qcond, выделяемая в вихревом проточном калориметре верхней охлаждаемой частью измерительной ТТ, оценивается следующим образом:
(ТсоМ - Тса1 )Fcond
Qcond 0
fcond + 0W li l w
1
1
(23)
Xl 1 w a cal acond
Значения параметров, указанных в скобках равны следующим значениям:
о fcond = 10-1 -10 -3мм; 0 w = 1мм;
acal = 2,4• 103 Вт/(м2 • K); acond > (5 -10)• 103 Вт /(м2 • K).
Введем общее выражение для теплового потока внутри ТТ вдоль продольной оси z:
дттт (z,tk)
dTTT (z ,tk)
zX-TT-—= z
д z
'TI
RTT ^TT (Z )
dz
Gvpr (Tb + GvpCvpd-df + K (1 - П ^^ (()Lev. (24)
dz
dz
dz
Решение уравнений (24) возможно только при положительных значениях времени и теплоемкости:
т>0; Сеу(Т)>0.
Для упрощения всех последующих расчетов можно принять начальное значение температуры равным нулю:
г ((,0 ) = 0 .
Уравнения (16) распространения теплового потока деу (г,1к) вдоль оси z в момент времени тк можно разделить и представить в виде системы из двух уравнений для расчета теплоемкости Сеук (г) и теплового потока деу(г,1к), значение которого связано с тепловым сопротивлением RTT и для момента времени тк может быть представлено следующим образом:
- м dqev (,tk .
zPevk ()t +--" = 0.
L
TT
dt (z,tk
RTT F {z
dz
d z
+ qev {Z,tk ) = 0
(25)
(26)
Значения переменной тепловой мощности деу (г ,1к) на внешней нижней плоской границе ТТ, на верхней границе плоского испарителя внутри парового канала и на верхней плоской границе ТТ равны:
qev { k ) = Qev; qev {0 ,035 ,tk ) = Q0,035 !
qev {-,tk) Qcond .
(27)
Уравнение (24) для определения передаваемой тепловой мощности капиллярно-пористым испарителем ТТ на высоте z=0,035ZТТ составлено с использованием распределения уменьшающейся температуры поверхности ТТ в адиабатическом калориметре:
qev [z = 0,035,tk ) = —Z-
= —0,035-
ZTT dt {z,t)
Rtt f {z) d z z=0,035
T2 LTT dt {z,t)
Rtt f {z) dz z=0,035
1. Проведение расчетов обратной задачи теплопроводности
Решение обратной задачи теплопроводности (ОЗТ) для испарителя ТТ представляет собой метод пошагового продолжения известного решения, например, равного нулю в начальный момент времени, для следующего временного шагаДтк = ( ,тк+1), соответствующего интервалу возрастания температуры при нагреве Дгк = (,гк+1) для расчета теплоемкости капиллярно-пористого испарителя высотой (толщиной) 0,035Трр.
Решение системы уравнений (25) и (26) осуществляется с использованием экспериментально определяемых значений температуры поверхности г (г ,Тк) в интервале области определения свойств испарителя ДТща0з5): (к к+1 )х(0;0,035 Ттт); экспериментальных значений теплового сопротивления Лрр (г ,Тк) в диапазоне всей длины ТТ ДТт: (,т к+1 )х(0,1Трт );и расчетных значений объемной теплоемкости испарителя Сеук (5?) также в интервале толщины испарителя (к ,т к+1) х (0 ;0,035 Трр) и удовлетворяющие уравнениям теплопроводности (16) и дополнительным условиям (19)-(22), при условии, что тепловое сопротивление, определенное в стационарных состояниях, известно как математический многочлен (14).
Расчетная схема Кранка-Николсона (см. рис. 5) с усреднением температуры на предыдущем и последующем временных расчетных слоях и с учетом температурной зависимости коэффициента теплопроводности может быть использована для повышения точности расчетов. Уравнение для расчета тепловой мощности испарителя по схеме Кранка — Николсона выглядит следующим образом:
peviCevi
ж Tfi1 — 2 Г
T.n+1 — Г-П
Dt
n+1 .
= 1e
тп+1 T+1
Ж rn+1
Tn
Ti—1 — 2Tin
(30)
T+1
2 h
2
2 h2
(28)
Коэффициент теплопроводности испарителя 1еи-был оценен ранее с помощью уравнения (31), все подробности расчетов приведены в работах [16]-[18]:
1 еу = П ■Х1 + (1 _ П)с. (31)
Уменьшающаяся с высотой тепловая мощность испарителя позволяет нам оценить работу ТТ. Для расчета теплоемкости испарителя Сегк в момент времени тк переменное выделение тепла в объеме испарителя можно рассматривать как внутренний процесс, интенсивность которого проявляет себя как меняющаяся температура внешней поверхности из уравнений (24) и (28). Минимальное значение тепловой мощности испарителя, которое было применено во всех расчетах, оценивают по формуле:
дде
0,035 Cevk {t )t = —
z ,tk
d z
(29)
z=0,035
На рис. 2 и 3 показано близкое к постоянному экспериментально определенное распределение температуры испарителя с кипящим диэтиловым эфиром при высокой тепловой нагрузке и дальнейшее более резкое снижение температуры в паровом канале ТТ.
Рис. 5. Шеститочечная разностная схема, в которой три точки берутся из следующего (нового) временного слоя n+1 и три из старого (предыдущего) временного слоя n Fig. 5. Six-point difference scheme in which three points are taken from the next (new) time layer n+1 and three ones —form the previous (old) time layer n
2
x
Экспериментальная температура tk в момент времени тк внутри парового канала ТТ с монотонным нагревом испарителя относительно температуры кипения tB диэ-тилового эфира (при атмосферном давлении) определяется следующим образом:
{к = ? (г ,тк ) = 5? (г ,тк) + ;
\ . (32)
% = 5?(г,тк) = (к - •
В начальный момент времени тк и при начальном значении температуры tk следует учитывать, что теплоемкость испарителя ТТ Сеук (5?) с начальным значением температурного напора 5?(г,тк) известна и для упрощения равна нулю, поэтому теплоемкость начиная с настоящего момента и при дальнейшем нагреве может быть представлена в виде полинома с переменной 5(0:
пс
Сеук () = Ф(Х)еСем (5?)'-1; пс < 10. (33)
1=1
Параметр аппроксимации Ф (X) является обязательной регуляризирующей величиной в виде монотонной возрастающей функции, Ф(Х)= X или, возможно, ехр(Х). Значения числового параметра Ф(Х) определяют стабильность и диапазон допустимых значений теплоемкости и сглаживают нелинейность при начале кипения в испарителе ТТ.
При переходе от линейной зависимости к экспоненциальной погрешность расчета теплоемкости заметно уменьшается, особенно когда полиномиальная функция теплоемкости выходит из области экстремума (области пика) (рис. 6). Это связано с тем, что ехр(Х) обеспечивает ограничение допустимых теплофизических свойств ТТ в положительном диапазоне значений: Сеук ,
,1"П7 > 0, и расчетное решение для Сеук значительно улучшается, особенно на выходе из пика теплоемкости, соответствующего фазовому переходу при начале кипения в испарителе ТТ.
В случае Ф(Х)= X погрешность восстановления теплоемкости в диапазоне температур 298^348 К (25^-75 °С). достигает 25%, а при замене на ехр(Х) падает до 2%.
Производная по времени от возрастающей температуры ТТ является постоянной величиной начиная с момента времени т > тк и является координирующим (синхронизирующим) параметром всех измерений в интервале времени Дтк:
¿1 > 0, I = 0...2 т>тк;
?к ( к )=?о+*к -3 -10 -3к/с;
Д?к = Дтк - 3 10 -3К/с.
,035 (к) = ?к0 ((= 0,035)0;
К ) = ?к0 ( = 0.
(34)
Экспериментальные значения температуры поверхности tk в момент времени тк вдоль образующей корпуса ТТ в вакуумной адиабатической оболочке известны, имеют положительный знак и могут быть записаны следующим образом:
?еу (Тк) = ?к0 ( = 0)
(35)
Учет и синхронизация равновесных значений теплового сопротивления Rттk (26) при медленном линейном
во времени нагреве и в моменты измерения времени тк выполняется с использованием значений увеличивающегося температурного напора 5tk из (32), (34), (35). Тепловая мощность из уравнения (24) рассчитывается следующим образом:
т>тк: 0=0 (тк )=^еу (0,тк);
0=0,035 (тк ) = Ъу (0,035,тк ); 0=1 (тк ) = 0сопё • (36)
Когда обратная задача (ОЗТ) в момент времени тк решена и температуры в момент времени тк+1 испарителя
г = 0:_ ?еу (%к) = ?ко ( = 0) на высоте г = 0,035: ,035 (тк) = = 0 ( = 0,035) и температура поверхности конденсации ТТ г = 1: ?1 (тк) = ?к0 ( = 1) измерены и рассчитаны, это решение (эксперимент и расчет) становится начальным распределением температуры и тепловой мощности при переходе к следующему временному интервалу Тк+1 ,ДТк+1 = (тк+1 ,тк+2). Из решения задачи на временном интервале Дтк+1 определяется функция Ск+1 (t) на температурном интервале 1к+1 = ((+1 ,?к+2). Чтобы уменьшить количество параметров Ск+1, определяемых на каждом временном интервале Дтк+1, перекрестную связь следует выполнять с помощью полинома (33) из предыдущего расчета в граничной точке температуры ^ по функции и ее производным, и это значительно упрощает все вычисления.
Наибольшее количество уравнений для сшивки N=2 (только теплоемкость или теплоемкость и ее производная по температуре) и максимальная степень полинома пс= 10 являются в первую очередь основными регуляризирую-щими параметрами расчета и были оценены на основе получаемой расчетной точности теплоемкости и экспериментальной погрешности измерения температуры.
Кроме того, в методе численного решения ОЗТ для коротких линейных ТТ с использованием полиномиального расширения теплоемкости, размер каждого временного шага Дтк также является важным регуляризующим параметром [12]—[15] и был выбран экспериментально. Характерное время опроса всех измерительных датчиков экспериментального стенда (технологический параметр), рис. 1, составляет Дтшеа!!=21 с, поэтому временной шаг Дтк при начале кипения диэтилового эфира в испарителе ТТ увеличился до Дтк=30 с, а температурный шаг уменьшился до значения Ы?к < 0,2 К. В начальной области нагрева температурный шаг может быть увеличен до значения Дтк <0,5 К при величине временного шага Дтк=60 с.
Программа, управляющая проведением измерений опытного стенда, (см. рис. 1), имеет два режима работы: «контроль» и «измерение». В режиме «контроль» циклически включаются все датчики, включая термопары, термисторы, емкостные датчики, контролируется мощность основного резистивного нагревателя и защитных нагревателей адиабатического калориметра, задается режим «нулевого нагрева» испарителя ТТ, результаты измерений обрабатываются и выводятся на экран дисплея.
В этом режиме задаются параметры управляющей программы, такие как длительность измерения показаний датчиков, время измерения цифровых вольтметров, осциллографа, частотомера и т. д. После того, как испа-
ритель ТТ и вихревой проточный калориметр достигают стационарных изотермических состояний, программа переключается в режим «измерение».
Главный резистивный нагреватель Н2 и система управления включаются и спустя 3 мин начинается линейный во времени нагрев испарителей ТТ, начинает работать адиабатическая система и все измерительные датчики. В этом режиме измеряется распределение температур на внешней поверхности измерительной ТТ, тепловая мощность испарителя, температура поверхности конденсации, тепловая мощность адиабатической системы и все другие тепловые характеристики, включая характеристики теплопередачи ТТ. Полученные массивы экспериментальных данных сохраняются и формируется библиотека измерений.
2. Расчет теплоемкости испарителя
Проанализируем случай, когда тепловое сопротивление RTT, уравнение (14) всей ТТ как целого известно, имеет положительное значение и необходимо рассчитать теплоемкость капиллярно-пористого испарителя высотой 3,5 мм, насыщенного диэтиловым эфиром.
Интегрируем уравнение (16) по г и, принимая во внимание экспериментальные граничные условия (19), (20), (21), мы получаем нелинейное интегральное уравнение для расчета теплоемкости испарителя в момент времени тк:
0,035 _ _ _ _
т гСеук (г)(г,тк)г = 0,035деУ (г,тк). (37)
0
Для решения таких интегральных уравнений используется метод конечных разностей, основанный на замене производных их приближенными значениями, выраженными в терминах конечных разностей вычисляемой функции в отдельных дискретных точках (узлах) конечно-разностной сетки, уравнение (30).
Решим функциональное уравнение (37) итерационным методом, для чего необходимо перейти к вариационной формулировке этой задачи. Мы определяем целевой функционал как расхождение, соответствующее разнице между левой и правой частями уравнения (37), и задача расчета для каждого временного интервала Дтк будет выглядеть как задача минимизации функционала расхождения:
!5Сеук () = 2е т
0,035 гРеук ()? (г ,тк )г-
т
0 _0,035 деу (к)
Л = шСук. (38)
Т =
Э15Сеук ( г ) к_
эСе
еук
0,035 Ж
т Л
гСеук (г)) +
0,035 ж
1 с
0
(г ,тк)+г
тт
Эг (г ^к)
RTTF (г) Эг
Эдеу (г ,тк)
эг
dгd т,
йг(1 %к +
(39)
где область переменных тк и деу (г,%к) сохраняется такой же, как внутри всего испарителя ТТ0 035 и всей длины ТТ;, а п и х являются множителями Лагранжа, соответствующими первому и второму уравнениям (25) и (26).
Чтобы получить градиентное выражение функционала L [32], вычислить множители Лагранжа и сформулировать сопряженную задачу расчета теплоемкости, мы представляем вариацию функции Лагранжа относительно переменных тк и деу (г ,Тк) и параметров теплоемкости С'еук (г), а также с учетом граничных условий (19) и (20) и начальных условий (25), (26), (34) и выражения для теплового потока (24) и (28):
тк+1
5Т = I Т
I т
0,035
Т %Сук () 5
тк+1
"1
тк
тк+1
тк
тк+1
С'еук ()Л + г Т
к
0,035 ж
1
0
0,035 X
Т (I - *
тт
Ятт F (г)
ц"
dzdТк_
тт
ЭЬ
5tdzdтк -
* (г) Эг
5деу (г ,тк )dzd тк + 1 глС'еук (+
(40)
0,035
1
- Т
0,035
тт
„ г./ 1_ 1 глС'еук Ог(,Тк_0,035деу (тк)
%т * ( ) 0
5td Тк.
Здесь величина gk обозначает производную от составляющих теплоемкости в следующем виде:
ЭСеук (г ,тк ).
' ЭСеук ((,тк)
Э15Сеук (г) тк+1
-= 1
ЭС'еук ((,тк) тк
0,035
| Лг(г,тк )
0
(41)
dгdтk,
где п — решение сопряженной задачи, полностью соответствующее исходной задаче (16), (19) и (21). При использовании обратного времени Дт=тк+! - тк — сопряженная задача записывается с помощью уравнений (25) и (26) и коэффициента теплопроводности ак:
Наиболее подходящим способом минимизации функционала является метод сопряженных градиентов [31]-[33]. Для вычисления градиентных составляющих функционала (38) (производных Фреше) в точках минимизации числовой последовательности используется метод теории оптимального управления процессом [32]-[33], описываемый уравнениями в частных производных [32]. В этом методе важную роль играют так называемые сопряженные краевые задачи, соответствующие исходной задаче оптимального управления. Функционал Лагранжа L для минимизации задачи (38) с условиями (16), (19), (20) и со значениями теплового потока на границах испарителя (25), (29) выглядит следующим образом:
гЛк = а Л— 1
Э(г-Элк / Эг
а =
Эг к
тт
Сеук ((,тк) Сеук ((,тк)ЯгТк*(г
(42)
Как начальные, так и граничные условия сопряженной задачи записывают обычным способом, и на верхней границе испарителя сопряженное значение равно калорическому остатку из квадратных скобок в уравнении (38):
л(( ,0 )= 0 ; л(0,Тк )=0.
т
+
к
т
+
к
+
г
к
2
т
к
к
к
1
h( = 0,035,tfc) =
,035 _
j ZCevk ((tk)d~Z-0,035qev (tk)
0,035 j
0
Z Lev■ 1 Qcond.
Э z
Таблица 1
Расчетные значения относительной теплоемкости при нагреве двух образцов с различными параметрами теплопередачи и оценка влияния различной математической обработки
Table 1
Design values of relative heat -absorbing capacity at heating two samples with different parameters of heat transfer and the estimation of different mathematical processing influence
Тип образца Ф (5) C/Co
Измерительный образец с теплопередачей теплопроводностью Ф (5)=5 3,3
Ф (5)=exp (5) 3,3
Модельный образец с теплопередачей за счет испарения жидкого диэтилового эфира и конденсации его пара Ф (5)=5 9,5
Ф (5)=exp (5) 11,8
_ э—
Z = 0 ■ ^"Э- = Qev - GeyQvP .
(45)
(tk+1 -t). (43)
Обратная задача для теплоемкости Сеук ((,тк) на каждом к -м интервале сводится к нахождению минимума функционала методом сопряженных градиентов.
Решение обратной задачи теплопроводности испарителя ТТ (25), (30) и (40) и решение интегрального уравнения (38) для определения теплоемкости испарителя было выполнено с использованием разработанной программы на Фортране [3]-[5]. В методе сопряженных градиентов одномерная минимизация в допустимом направлении осуществлялась с использованием метода интерполяции-экстраполяции [32]. Описанный алгоритм определения кусочной функции Сеук(5?) с известной функцией теплового сопротивления $тт
(?) ранее был
протестирован на нескольких классических моделях обратных задач теплопроводности. Полученные результаты решения этих задач показывают, что функция теплоемкости Сеук (5?) может быть восстановлена во всех случаях с погрешностью менее 5%.
В табл. 1 приведены результаты решения обратной задачи теплопроводности для двух цилиндрических образцов, а именно измерительного образца с теплопередачей теплопроводностью и модельного образца с локальным фазовым переходом, с одинаковой длиной образов 3-10-3 м. Теплоемкость составляет С0 = 1,106 Дж/(м3-К), тепловое сопротивление Катр1е = 0,25 -10 К/Вт , теплота испарения диэтилового эфира 0ур = 345 кДж/кг при температуре испарителя Тгу=328,75 К (55,6 °С) и давлении Р~6,13-105 Па , 0у = 210 Вт , скорость нагрева эфира 3 -10-3 К/с.
Для измерительного образца с теплопередачей теплопроводностью все подробности взяты из уравнения (16). Для модельного образца граничные условия приведены ниже в (44) и (45):
Здесь скорость испарения и давление насыщенного пара оцениваются следующим образом:
t ( - P) Gev, = const—7= ' ;
V2 pR—M
Pev= P • exp
Qvp_ R—
Вся численная реализация представленного алгоритма была выполнена путем замены всех этих уравнений их разностными аналогами, интегралы были рассчитаны с использованием известного метода Симпсона [34]-[36] и численной схемы Кранка — Николсона, упомянутой выше в (30) и показанной на рис. 5. Расчетная схема Кранка-Николсона требует примерно в 10 раз меньше вычислительной работы для ее реализации при расчете на одном и том же временном интервале [34]-[35].
Итерационный процесс завершают при получении результата 5Сеу/Сеу, не превышающего случайной погрешности скорости роста температуры и тепловой мощности:
5Cev= З dt f +
C,
d9ev Чеу
(46)
Результаты расчета и решения уравнения (37) в диапазоне температурного напора на испаритель ТТ 5? = Теу - ТБ = (-5... + 25) К показаны на рис. 6, случайная погрешность при расчете теплоемкости не превышает (2-3) %.
На графике: Сеу — теплоемкость насыщенного ди-этиловым эфиром испарителя, Дж/К; Сеу0 — теплоемкость испарителя опорной ТТ, заполненного осушенным воздухом, Дж/К; 1 — значения относительной теплоемкости нижнего фрагмента короткой ТТ (фрагмента испарителя) с выполненным в виде сопла Лаваля паровым каналом,
(44)
Рис. 6. Расчетное значение теплоемкости испарителя ТТ dCev/Cev0
Fig. 6. Deign value of heat-absorbing capacity for ТТ evaporator dCevlCev0
2
2
t
полученные путем решения обратной задачи теплопроводности с шагом по температуре 0,5 К 5/= Теу -Тв, при начинающемся процессе кипения диэтилового эфира; 2 — полином десятой степени для сглаживания полученных точек теплоемкости испарителя.
Описанный численный алгоритм решения обратной задачи теплопроводности коротких ТТ и расчета теплоемкости испарителя ТТ с использованием экспериментально определенных значений распределения температуры по поверхности внутри адиабатической оболочки и теплового сопротивления Лтт (5г) в зависимости от температурной нагрузки на испаритель 5г = Теу _ Тв и тепловой мощности позволяет оценить параметры и характеристики кипящего испарителя [37]-[39]. Вычисленное значение теплоемкости можно проверить следующим образом. Энтальпия испарения работающего испарителя с кипящим диэтиловым эфиром может быть оценена по уравнению (47), в которое подставлены известные табличные данные для диэтилового эфира и получено выражение для эффективной теплоемкости нагреваемого капиллярно-пористого испарителя ТТ при начале кипения эфира:
С = Ну _(345 +1,75 -11 + 2,34-16)1,57 -10_4) =
Сеу : С -
г Сеу0
= 21,02 кДж/К (47)
Отношение теплоемкости испарителя с диэтиловым эфиром к теплоемкости сухого испарителя 5Сеу/Сеу0 с учетом теплоемкости резистивного нагревателя Сн2 = 1,2 кДж/К выглядит следующим образом, что очень близко к максимальному численному значению теплоемкости испарителя ТТ:
Сеу _Сеу0 + Неу / г Ну / / = 21,02 + 1,2
С
еу0
С
еу0
С
еу0
С
еу0
=18,5. (48)
Уравнения (47) и (48) подтверждают соответствие экспериментальных и расчетных теплофизических характеристик диэтилового эфира, находящегося в испарителе при высоком температурном напоре и кипении эфира.
Выводы
1. Представленный расчетный результат экстремального поведения теплоемкости испарителя коротких ТТ с многослойной капиллярно-пористой вставкой (фитилем), изготовленной из металлической сетки с различными размерами ячеек, образующей внутренний паровой канал в виде сопла, близкого к соплу Лаваля, с большой подачей диэтилового эфира из пор вставки в плоский капиллярно-пористый испаритель с инжекторными каналами, при высоком температурном напоре на испаритель 5г = Теу _ Тв > 11 К, полученный путем решения обратной задачи теплопроводности (ОЗТ) с использованием экспериментально определенного распределения температуры поверхности корпуса ТТ вдоль образующей и тепловой мощности Qeу 035,0со^ при монотонном и линейном во времени нагреве, подтверждает возникновение процесса кипения в испарителе.
2. Температура внешней поверхности стенки ТТ вокруг капиллярно-пористого испарителя, близкая к постоянному значению при монотонном нагреве, также
подтверждает начало процесса кипения и позволяет оценить комплекс теплофизических параметров диэтилово-го эфира в испарителе.
3. С помощью пошаговой регуляризации решения обратной задачи теплопроводности ТТ и кусочно-гладкой полиномиальной аппроксимации неизвестных коэффициентов удается уменьшить ошибку численного восстановления теплоемкости испарителя, используя функцию полиномиально-экспоненциальной аппроксимации, которая ограничивает диапазон допустимых значений теплоемкости только положительными значениями.
Обозначения, применяемые в статье
а — коэффициент температуропроводности, м2/с; ак — к-е значение температуропроводности ТТ, м2/с; СТТ (Г) — теплоемкость ТТ как целого, Дж/К; е„ (Г) — объемная теплоемкость ТТ, Дж/ (м3 К); е„ (Т) — удельная изохорная теплоемкость ТТ, Дж/ (кгК); Сеу0 — теплоемкость сухого капиллярно-пористого испарителя, Дж/К;
Сеу — теплоемкость насыщенного диэтиловым эфиром капиллярно-пористого испарителя, Дж/К; Сеу — вычисленная теплоемкость испарителя, насыщенного
диэтиловым эфиром, Дж/К; С — теплоемкость диэтилового эфира, Дж/ (кгК); Ср — теплоемкость влажного пара диэтилового эфира, Дж/ (кгК);
Сйу- теплоемкость жидкого диэтилового эфира на поверхности
испарителя, Дж/ (кгК); ЕН2 — тепловая мощность резистивного нагревателя, Вт; ЕНР (2) — площадь поперечного сечения ТТ, м2; Gl — массовый поток жидкой фазы диэтилового эфира в испаритель ТТ, кг/с; Gшix — массовый поток влажного пара диэтилового эфира над
испарителем ТТ, кг/с; Оур — массовый поток сухого пара диэтилового эфира над испарителем ТТ, кг/с; (1 -xev)=G'/ (G'+G) — степень влажности влажного пара в ТТ; gk — производная компонент теплоемкости; Яеу — энтальпия испарителя ТТ, Вт; ЬГТ — длина ТТ, м;
Ьеу—толщина капиллярно-пористого испарителя, м; Р — среднее давление в паровом канале ТТ, Па; Реу — давление над испарителем в ТТ, Па; РсС0пЛ — давление на поверхности конденсации, Па; Ру — давление внутри парового тороидального вихревого кольца, близкое к давлению конденсации Р^^; Рг — число Прандтля пара диэтилового эфира; Qev — тепловая мощность, поступающая в испаритель ТТ, Вт; О.с0пс1 — тепловая мощность конденсации, W; Rc — тепловое сопротивление слоя смазки на нижнем торце
испарителя, К/Вт; Rev — тепловое сопротивление испарителя, К/Вт; Rw — тепловое сопротивление металлических крышек, К/Вт; RТТ (Т) — тепловое сопротивление ТТ целиком, К/Вт; г(Тв) — удельная теплота испарения диэтилового эфира, кДж/кг; Тв — температура кипения диэтилового эфира, К; Тса1 — температура охлаждающей воды в проточном калориметре, К;
Тс0п1) — температура поверхности конденсации ТТ, К;
T fcond — средняя температура слоя диэтилового эфира на поверхности конденсации ТТ, К; T fev — средняя температура кипящего слоя диэтилового эфира на поверхности испарителя, К; rTT(z) — температура внутри парового канала ТТ на высоте
z над испарителем, K; Tev — стационарная температура сетчатого испарителя, насыщенного диэтиловым эфиром, K; Tsc — стационарная температура металлической сеточной
структуры испарителя, K; Tv — температура внутри тороидального парового вихря, K; t — скорость роста температуры испарителя, K/c; v — скорость движения паровой и жидкой фаз внутри испарителя, м/с;
xev — степень влажности пара в испарителе; - = z / LHP — безразмерная продольная координата внутри ТТ; uv, vv — аксиальная и тангенциальная компоненты скорости
вихря в паровом канале ТТ, м/с; acal — коэффициент теплоотдачи в вихревом проточном калориметре, Вт/ (м2К); acond — коэффициент теплоотдачи в области конденсации ТТ, Вт/ (м2-К);
8Jiond — средняя тощина пленки диэтилового эфира на поверхности конденсации, м; Sw — толщина нижней металлической крышки ТТ, м; П — множитель Лагранжа;
X j — коэффициент эффективной теплопроводности ТТ, Вт/ (м-К);
l HP — коэффициент средней теплопроводности ТТ, Вт/ (мК); Xj — коэффициент теплопроводности диэтилового эфира, Вт/ (мК);
X (t) — коэффициент теплопроводности ТТ, Вт/ (мК); Xsc — коэффициент теплопроводности металлического сеточного испарителя, Вт/ (мК); Xw — коэффициент теплопроводности верхней и нижней металлических крышек ТТ, Вт/ (м-К); ^ — коэффициент трения; £, — коэффициент конденсации ТТ; П — пористость испарителя и вставки ТТ; РТТ — средняя плотность ТТ, кг/м3; Р7 — плотность влажного пара внутри ТТ, кг/м3; pv — плотность пара внутри тороидального парового вихря, кг/м3;
-ev — средняя плотность испарителя, кг/м3; т — время, с.
Дтк — временной шаг при расчетах и измерениях, с; Ф^ — диссипативная функция; X — множитель Лагранжа.
Литература/References
1. Brebbia C. A., Domínguez J. Boundary Elements. An Introductory Course. WIT Press. 1998. 325 p.
2. Nakayama A., Kuwanara F., A general bioheat transfer model based on the theory of porous Media. Int. J. of Heat and Mass Transfer. 2008. v. 51, pp 3190-3199.
3. Beck J., Blackwell B., Clair C. Inverse heat conduction 111 posed problems. NY. 1985. 308 p.
4. Карслоу Г., Егер Д. Теплопроводность твердых тел. М.: Наука 1964. 488 с.
5. Лыков А. В. Теория теплопроводности М.: Высшая школа, 1967, 600 с. [Lykov A. V. Theory of thermal conductivity Moscow, Higher School, 1967, 600 p. (in Russian)]
6. Платунов Е. С. Теплофизические измерения в монотонном режиме. 1973. Л.: Энергия, 144 с. [Platunov E. S. Thermo-physical measurements in monotone mode. 1973. Leningrad: Energiya, 144 p. (in Russian)]
7. Исаченко В. П., Осипова В. А., Сукомел А. С. Теплопередача. М.: Энергоиздат, 1981. 416 с. [Isachenko V. P., Osipova V. A., Sukomel A. S. Heat transfer. Moscow: Energoizdat, 1981. 416 p. (in Russian)]
8. Кутателадзе С. С. Основы теории теплообмена. М.: Ато-миздат, 5-е изд. 1979. 416 с. [Kutateladze S. S. Fundamentals of the theory of heat transfer. Moscow, Atomizdat, 5th ed. 1979. 416 p.]
9. Денисов А. М. Введение в теорию обратных задач. М.: МГУ. 1994. 208 с. [Denisov A. M. Introduction to the theory of inverse problems. Moscow: MSU. 1994. 208 p. (in Russian)]
10. Тихонов А. М. О решении некорректно поставленных задач и методе регуляризации // ДАН СССР 1963. т. 151. № 3. с. 501-504. [Tikhonov A. M. On solving incorrectly set tasks and the regularization method // DAN USSR 1963. vol. 151. No. 3. pp. 501-504. (in Russian)]
11. Тихонов А. Н., Арсенин В. Я. Методы решения некорректных задач. М.: Наука, 1979. 284 с. [Tikhonov A. N., Arsenin V. Ya. Methods of solving ill-posed problems. Moscow: Nauka, 1979. 284 p. (in Russian)]
12. He J.-H. Variational iteration method — a kind of non-linear analytical technique: some examples. International Journal of Non-Linear Mechanics, 1999, v. 34, N. 4, pp. 699-708.
13. He J.-H. Some asymptotic methods for strongly nonlinear equations. International Journal of Modern Physics B, 2006. v. 20, N. 10, pp. 1141-1199.
14. Lu J. Variational iteration method for solving two-point boundary value problems. Journal of Computational and Applied Mathematics, 2007. v. 207, N. 1, pp. 92-95.
15. Yan L., Fu C.-L., Yang F. L. The method of fundamental solutions for the inverse heat source Problem. Engineering Analysis with Boundary Elements, 2008. v. 32, N. 3, pp. 216-222.
16. Seryakov A. V. Intensification of heat transfer processes in the low temperature short heat pipes with Laval nozzle formed vapour channel. American Journal of Modern Physics. 2018. v. 7, N 1, pp. 48-61.
17. Seryakov A. V, Shakshin S. L., Alekseev And. P. The droplets condensate centering in the vapour channel of short low temperature heat pipes at high heat loads. Journal of Physics: Conference Series. 2017. 891. 012123.
18. Seryakov A. V. Computer modeling of the vapour vortex orientation changes in the short low temperature heat pipes. I. J. Heat and Mass Transfer. 2019. v. 140. pp. 243-259.
19. Серяков А. В. Исследование характеристик коротких низкотемпературных тепловых труб с паровым каналом в виде сопла. // ПМТФ. 2016. Т. 57, № 1, с. 80-94. [Seryakov A. V. Investigation of the characteristics of short low-temperature heat pipes with a steam channel in the form of a nozzle. PMTF. 2016. vol. 57, No. 1, pp. 80-94. (in Russian)]
20. Seryakov A. V. The study of condensation processes in the low-temperature short heat pipes with a nozzle-shaped vapour channel. Engineering 2017. v. 9, pp. 190-240.
21. Seryakov A. V. Numerical modeling of the vapour vortex. J. high energy physics, gravitation. 2019, N 5, pp. 218-234.
22. Fowler, R. H., Statistical Mechanics, 2nd Ed., London, 1936. Cambridge University Press.
23. Chapman, S. and Cowling, T. G. The Mathematical theory of non-uniform gases. London, Cambridge University Press. 447 p. 1953.
24. CFdesign 10.0 2009. Version 10.0-20090623. User's Guide.
25. Handbook of Physical-Chemical Properties and Environmental Fate for Organic Chemicals. Vol. III. Oxygen containing compounds. D. Mackay, W. Shiu, K. Ma, S. Lee. Taylor and Francis. 2006.
26. Таблицы физических величин: Справочник. Под ред. Кикоина. М.: Атомиздат 1976. 1008 с. (in Russian)]
27. Варгафтик Н. Б. Справочник по теплофизическим свойствам газов и жидкостей. М.: Наука, 1972. 721 с. [Varga-ftik N. B. Handbook of thermophysical properties of gases and liquids. Moscow: Nauka,1972. 721 p. (in Russian)]
28. Faghri A. Heat Pipe Science and Technology. Washington USA, Taylor and Francis. 1995.
29. Odibat Z. M. and S. Momani. Application of variational iteration method to nonlinear differential equations of fractional order. International Journal of Nonlinear Sciences and Numerical Simulation, 2006. v. 7, N. 1, pp. 27-34.
30. Farcas A., Lesnic D. The boundary-element method for the determination of a heat source dependent on one variable. Journal of Engineering Mathematics, 2006. v. 54, N. 4, pp. 37388.
31. Эльсгольц Л. Э. Вариационное исчисление: учебник для физических и физико-математических факультетов университетов. 7-е изд. М.: 2008. 205 с. [Elsholts L. E. Calculus of variations: textbook for physical and physical-mathematical faculties of universities. 7th ed. Moscow: 2008. 205 p. (in Russian)]
32. Эльсгольц Л. Э. Дифференциальные уравнения и вариационное исчисление. 2-е изд. М.: Наука, 1969. 424 c.
[Elsholts L. E. Differential equations and calculus of variations. 2nd ed. Moscow: Nauka, 1969. 424 p. (in Russian)]
33. Цлаф Л. Я. Вариационное исчисление и интегральные уравнения: справочное руководство. 3-е изд. Санкт-Петербург. 2005. 192 с. [Tslaf L. Ya. Calculus of variations and integral equations: a reference guide. 3rd ed. St. Petersburg. 2005. 192 p. (in Russian)]
34. Марчук Г. И. Методы вычислительной математики. М.: Наука. Главная редакция физико-математической литературы, 1980. 536 с. [Marchuk G. I. Methods of computational mathematics. M.: Nauka. The main edition of the physical and mathematical literature, 1980. 536 p. (in Russian)]
35. Тихонов А. Н., Самарский А. А. Уравнения математической физики. М.: Наука, 1966. 724 с. [Tikhonov A. N., Samar-skiy A. A. Equations of mathematical physics. M.: Nauka, 1966. 724 p. (in Russian)]
36. Калиткин Н. Н. Численные методы. М.: Наука, 1978. 512 с. [Kalitkin N. N. Numerical methods. M.: Nauka, 1978. 512 p. (in Russian)]
37. Seryakov A. V., Alekseev A. P. A study of the short heat pipes by the monotonic heating method Journal of Physics: Conference Series. 2020. 1683. 022051.
38. Seryakov A. V., Alekseev A. P. Investigation of short heat pipes by the method of monotonous Heating. III International Conference «Modern problems of Thermophysics and Energy», Moscow, NRU MEI, 2020 October 19-23, pp. 270-272.
39. Серяков А. В., Конькин А. В., Алексеев А. П. Экспериментальное исследование тепловых характеристик линейных тепловых труб при воздействии внешних вибраций // Вестник Международной академии холода. 2020. № 3. С. 27-38. [Seryakov A. V., Konkin A. V., Alekseev A. P. Experimental study of linear heat pipes thermal characteristics under the external vibrations effect. Journal of International Academy of Refrigeration. 2020. No 3. p. 27-38. (in Russian)]
Сведения об авторах
Серяков Аркадий Владимирович
К. ф.-м. н., ООО «Рудетранссервис», 173003,
г Великий Новгород, ул. Большая Санкт-Петербургская, 64,
seryakovav@yandex.ru
Алексеев Андрей Петрович
Член совета директоров, ООО «Рудетранссервис», 173003, г Великий Новгород, ул. Большая Санкт-Петербургская, 64
Information about authors
Seryakov Arkady V.
Ph. D., LLC «Rudetransservice», Sankt-Petersburg Str. 64, Veliky Novgorod, 173003, Russia, seryakovav@yandex.ru
Alekseev Andrey P.
BoD member, LLC «Rudetransservice», Sankt-Petersburg Str. 64, Veliky Novgorod, 173003, Russia
