Теплофизика высоких температур, 2020, T. 58, № 1, стр. 113-122

ВВЕДЕНИЕ

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

Конкретная прикладная задача может приводить к дифференциальному уравнению любого порядка или к системе уравнений любого порядка. Однако известно, что обыкновенное дифференциальное уравнение N-го порядка [1]

при помощи замены z^(k)(τ) ≡ z_k(x) можно свести к эквивалентной системе N уравнений первого порядка [1] где z₀(τ) = z(τ). Аналогично произвольную систему дифференциальных уравнений любого порядка можно заменить некоторой эквивалентной системой уравнений первого порядка. Поэтому в дальнейшем будет рассматриваться система обыкновенных дифференциальных уравнений первого порядка [1] или в векторной форме

Для определения вектора z достаточно начального условия (задача Коши)

Такими математическими моделями с сосредоточенными параметрами (МСП), основанными на обыкновенных дифференциальных уравнениях, и описываются теплофизические системы [2]:

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

Выражая из (2) T_i, _{j+ 1}, получаем квазилинейную нестационарную дискретную динамическую систему следующего вида:

где x_i = (T_{i, 0}, T_{i, 1}, …, T_i,M) – вектор состояния, f = = (f₁, f₂, …, f_N).

Рассмотрен детерминированный подход к параметрической идентификации такой системы, определены условия разрешимости поставленной задачи, сформулирован критерий существования решения, и построен алгоритм его поиска, возможности которого продемонстрированы на численных примерах [4, 5].

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

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

В [9] исследовано влияние внутреннего тепловыделения и числа Био на точность и применимость моделей с сосредоточенными параметрами в нестационарных одномерных задачах теплообмена. Для получения средней температуры и плотности теплового потока в пространственном направлении применены классическая МСП на основе объемной теплоемкости и усовершенствованная модель на основе эрмитовой и полиномиальной аппроксимации. Введены количественные критерии оценки применимости моделей с сосредоточенными параметрами.

В [10] технологическая система сложного вида рассмотрена как “черный ящик” со входом и выходом. Описан тепломассоперенос в этой системе при помощи обыкновенных дифференциальных уравнений с постоянными коэффициентами определенного порядка. Предложена методика использования “неустойчивых” уравнений для описания практических задач. Под неустойчивостью понимается тот факт, что полученный дифференциальный оператор имеет характеристические числа с положительной вещественной частью. Определены коэффициенты уравнения на основе анализа кривой отклика системы на возмущение в виде “ступеньки”. Рассмотрен конкретный пример реализации данного алгоритма.

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

Метод итерационной регуляризации адаптирован для математических моделей теплообмена с сосредоточенными параметрами: построены аналитические выражения для расчета неопределенных множителей Лагранжа при определении градиента функционала невязки температуры и расчета шага спуска, работоспособность которых подтверждена численными экспериментами. Полученные результаты обобщены для технических систем космических аппаратов [12, 13].

Анализ вариационного метода регуляризации [14] и метода итерационной регуляризации [15–17] позволил разработать модифицированный вариационный метод итерационной регуляризации (итерационно-вариационный метод) решения задач идентификации математических моделей с сосредоточенными параметрами на основе минимизации сглаживающего функционала и функционала температурной невязки [18, 19]. Полученные результаты применены для идентификации математических моделей составных частей космического аппарата (КА), в частности теплового сопротивления и излучательной способности экранно-вакуумной тепловой изоляции (ЭВТИ) [19, 20].

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

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

ЭКВИВАЛЕНТНАЯ ТЕПЛОВАЯ СХЕМА И МАТЕМАТИЧЕСКАЯ МОДЕЛЬ С СОСРЕДОТОЧЕННЫМИ ПАРАМЕТРАМИ СОТОВОЙ ПАНЕЛИ И ЭКРАННО-ВАКУУМНОЙ ТЕПЛОИЗОЛЯЦИИ

Рассмотрим характерный элемент конструкции КА в виде многослойной теплопередающей стенки, состоящей из последовательно расположенных ЭВТИ и СП. СП устанавливается на силовом каркасе и вместе с ним образует несущую конструкцию КА.

Расчет проводимостей ЭВТИ и СП проводится по показаниям двух групп аналоговых датчиков температуры (АДТ) в режимах тепловых вакуумных испытаний (ТВИ), имитирующих предельные условия эксплуатации (рис. 1): “Д1” – внешний нагрев инфракрасными имитаторами (ИКИ); “Д2”, “Д4” – охлаждение криоэкранами с температурой 83 К; “Д3” – внутренний нагрев опытного образца за счет тепловыделений аппаратуры и электронагревателей.

Прямой расчет проводимостей невозможен из-за отсутствия экспериментальных значений плотности теплового потока, проходящего через многослойную стенку, поэтому проводимости ЭВТИ и СП уточняются решением задачи идентификации МСП в удельных величинах, включающих три узла, соответствующих физическим поверхностям: узел 7 – внешняя поверхность ЭВТИ, узел 8 – внешняя поверхность СП, узел 9 – внутренняя поверхность СП (рис. 2).

Узел 7 – модельная оболочка ЭВТИ толщиной 2.59 мм (толщина ЭВТИ в сжатом состоянии), имеющая с изолируемой поверхностью тепловой контакт, проводимость которого соответствует проводимости ЭВТИ в поперечном направлении. Плотность модельной оболочки равна

где m_ЭВТИ = 0.4 кг – масса ЭВТИ площадью F = 1 м²; δ_ЭВТИ = 2.59 × 10^–3 м – толщина ЭВТИ в сжатом состоянии. Удельная теплоемкость оболочки ЭВТИ принята равной с₇ = с_ЭВТИ = 200 Дж/(кг К).

Узлы 8, 9 моделируют две алюминиевые пластины толщиной

между которыми задана тепловая проводимость, эквивалентная теплопроводности сот в поперечном направлении. Плотность пластин пересчитывается с учетом их массы:

Удельная теплоемкость модельных поверхностей СП равна

Удельное сопротивление ЭВТИ имеет начальное значение

а удельное сопротивление СП где разности температур определены экспериментально.

Соответствующие удельные проводимости:

Система уравнений математической модели с сосредоточенными параметрами записывается в следующем виде:

где q_{s, 7} = ε₇q₇ – плотность поглощенного поверхностью ЭВТИ инфракрасного потока q₇, Вт/м²; q₉ – плотность теплового потока, поступающего в СП от внутренних источников тепла, Вт/м².

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

Будем искать излучательную способность ЭВТИ ε₇ = ε₇(τ) и удельные тепловые проводимости ЭВТИ α₇₈ = α₇₈(τ) и СП α₈₉ = α₈₉(τ) как функции времени на временнóй сетке τ_j (j = 0, 1, 2, …, M + 1): ε_j = = ε₇(τ_j), α_78j = α₇₈(τ_j), α_89j = α₈₉(τ_j). Тогда система уравнений (3) запишется следующим образом:

с начальным условием и условием непрерывности

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

Производные функций равны

Тогда выражения метода Ньютона и соответствующие им приращения искомых функций запишутся как

где l = 1, 2, …; $\varepsilon _{7}^{{(0)}}$ = 0.01, $\alpha _{{78}}^{{(0)}}$ = 0.1 Вт/(м² К), $\alpha _{{89}}^{{(0)}}$ = = 1.565 Вт/(м² К) – начальные приближения.

Для регуляризации решения используется дополнительная информация в виде экспериментальных температур

Регуляризированное решение выражается через решение на текущей итерации l, его приращение относительно итерации (l – 1) и безразмерный шаг спуска β_ij (i = 7, 8, 9; j = 1, 2, …, M):

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

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

Тогда температуры T_ij изменятся на некоторую величину υ_ij(τ) (i = 1, …, N;  j = 1, …, M), а вариация температуры υ_ij(τ) удовлетворяет краевой задаче:

с начальным условием и условием непрерывности

Функционал температурной невязки на итерации l запишется как

Следовательно, приращения компонент функционала равны

а вариации температур υ_ij определяются из следующего выражения:

Подставляя вариации температуры из (9) в (7), получаем уравнение для определения β_ij (i = 7, 8, 9; j = 1, 2, …, M):

где температуры T_ij (i = 7, 8, 9; j = 1, 2, …, M) являются решением краевой задачи (4)–(6) на текущей итерации, а производные dT_ij/dτ выражаются из системы уравнений (4) после ее разрешения относительно температур T_ij. Выражая производные dT_ij/dτ из уравнений (4), подставляя их в уравнения (10) и решая их относительно β_ij, получаем следующие выражения безразмерного шага спуска:

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

ЭКСПЕРИМЕНТАЛЬНЫЕ ДАННЫЕ И РЕЗУЛЬТАТЫ ИДЕНТИФИКАЦИИ

Экспериментальные температуры узлов математической модели показаны на рис. 3, плотность потока излучения от инфракрасных имитаторов – на рис. 4, а плотность теплового потока, имитирующего внутреннюю тепловую мощность в узле 9, – на рис. 5.

Дисперсия экспериментальной температуры $T_{i}^{{(0)}}$ равна

где φ_ij = $T_{{ij}}^{{\left( 0 \right)}}$ – $T_{{ij}}^{'},$ $T_{{ij}}^{'}$ – наиболее близкое к истинному значение температуры, чаще всего среднее арифметическое ряда измерений в этой точке. Так как данный режим испытаний является единственным, определить среднее арифметическое температуры в данный момент времени невозможно, поэтому приближенно можно считать φ_ij ≈ $\Delta T_{{ij}}^{{\left( 0 \right)}},$ где $\Delta T_{{ij}}^{{\left( 0 \right)}}$ = 2 К – абсолютная погрешность измерения температуры по показаниям датчика. Тогда

Следовательно, дисперсия температуры одинакова для всех трех узлов математической модели, но изменяется в зависимости от количества M = (M₁, M₂, M₃, M₄) = (25, 37, 69, 31) точек временнóй сетки в различных режимах: D_i = (D_i1, D_i2, D_i3, D_i4) = (0.167, 0.111, 0.059, 0.133) К² (i = 7, 8, 9; второй индекс определяет номер режима). Среднеквадратичная ошибка температурных измерений для всей системы равна

т.е. в зависимости от экспериментального режима δ_i = (δ_i1, δ_i2, δ_i3, δ_i4) = (0.501, 0.333, 0.177, 0.399) К².

Для проверки устойчивости решения на распределения экспериментальных температур накладывается нормальное возмущение со средним квадратичным отклонением 0.4 К, соответствующим дисперсии температуры 0.16 К². Устойчивость решений проверяется сопоставлением функций, определенных по возмущенному температурному полю, с функциями, идентифицированными по исходным распределениям экспериментальных температур.

Распределение значений температурного функционала (8) в зависимости от номера итерации показано на рис. 6. Наименьшие значения функционала, соответствующие наилучшему решению, достигаются (в скобках для возмущенного решения): в режиме “Д1” – 763.2 К², на первой итерации (775.3 К², на первой итерации); в режиме “Д2” – 376.8 К², на четвертой итерации (418.1 К², на четвертой итерации); в режиме “Д3” наименьшее значение функционала, равное 676.0 К², достигается на двенадцатой итерации, но после девятой итерации со значением функционала 1362 К² гладкость решения ухудшается, поэтому решение считается оптимальным на девятой итерации (4282 К², на девятой итерации); в режиме “Д4” – 374.0, на восьмой итерации (454.0 К², на тринадцатой итерации).

Начальные приближения задавались следующие (в скобках для возмущенного решения): в режимах “Д1”, “Д2” – $\varepsilon _{7}^{{\left( 0 \right)}}$ = 0.01, $\alpha _{{78}}^{{\left( 0 \right)}}$ = 0.1 Вт/(м² К), $\alpha _{{89}}^{{\left( 0 \right)}}$ = 1.0 Вт/(м² К) (те же); в режиме “Д3” – $\varepsilon _{7}^{{\left( 0 \right)}}$ = = 0.008815 (0.009587), $\alpha _{{78}}^{{\left( 0 \right)}}$ = 0.03866 (0.04246) Вт/(м² К), $\alpha _{{89}}^{{\left( 0 \right)}}$ = 0.1452 (0.07836) Вт/(м² К); в режиме “Д4” – $\varepsilon _{7}^{{\left( 0 \right)}}$ = 0.01268 (0.01187), $\alpha _{{78}}^{{\left( 0 \right)}}$ = 0.05503 (0.04915) Вт/(м² К), $\alpha _{{89}}^{{\left( 0 \right)}}$ = 0.3948 (0.2002) Вт/(м² К). В режимах “Д1”, “Д2” начальные приближения определены по предварительной теоретической оценке, в режимах “Д3”, “Д4” – по последней временнóй точке предыдущего участка для обеспечения лучшей гладкости решения на границе режимов.

Идентифицированные функции представлены на следующих рисунках: излучательная способность ЭВТИ ε₇ – на рис. 7, удельное тепловое сопротивление ЭВТИ R₇₈ – на рис. 8, удельное тепловое сопротивление СП R₈₉ – на рис. 9.

Излучательная способность ЭВТИ имеет глобальный минимум 0.00446 (0.00458) в режиме “Д3” и глобальный максимум 0.0389 (0.03883) в режиме “Д4”, на два–три порядка отличаясь от излучательной способности наружного слоя (0.63) и являясь некоторой суперпозицией свойств внутренних слоев. Значение ε₇(τ) зависит не только от температуры T₇ наружного слоя, но и от температуры T₈ внутреннего слоя, а также от разности этих температур и теплового сопротивления R₇₈, т.е. от величины и направления плотности теплового потока. В результате, характер излучательной способности ЭВТИ сложен и нелинеен.

Тепловое сопротивление ЭВТИ R₇₈ имеет глобальный минимум 6.2 (6.2) м² К/Вт в режиме “Д4” и глобальный максимум 55.6 (54.3) м² К/Вт в режиме “Д3”. Глобальные и локальные максимумы теплового сопротивления ЭВТИ соответствуют глобальным и локальным минимумам излучательной способности ЭВТИ и наоборот.

Тепловое сопротивление СП R₈₉ имеет глобальный минимум 0.411 (0.397) м² К/Вт в режиме “Д2” и глобальный максимум 9.4 м² К/Вт в режиме “Д4” (12.8 м² К/Вт в режиме “Д2”). Качественно поведение теплового сопротивления СП совпадает с характером временнóй зависимости термического сопротивления ЭВТИ, имея при этом точку разрыва на границе режимов “Д2”, “Д3”. По всей видимости, наличие точки разрыва объясняется изменением величины и направления теплового потока через сотовую панель в результате включения внутренних источников тепла. Причем если теоретическое значение плотности теплового потока от внутренних источников изменяется по экспоненте (рис. 5), то экспериментальная тепловая мощность практически сразу же достигает номинального уровня.

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

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

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

Итерационный процесс продолжается до достижения минимального значения температурного функционала. Чаще всего ему соответствует наиболее гладкое решение. Тем не менее возможна ситуация, когда минимуму функционала соотносится решение с неприемлемой дисперсией, как, например, в режиме “Д3”. Тогда необходимо найти итерацию с минимальной дисперсией решения, после которой она возрастает, и определить наилучшее соотношение между гладкостью решения (значением дисперсии) и значением температурного функционала.

В итерационной регуляризации не всегда удается достичь значения функционала, сопоставимого по точности с исходными данными. Итерационный процесс останавливается по условиям достижения минимума функционала или замедления его сходимости к минимуму. В редких случаях (функционал считается квазимонотонным) возможно получение нескольких локальных минимумов и сравнение их между собой.

ЗАКЛЮЧЕНИЕ

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

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

3. Излучательная способность ЭВТИ изменяется от 0.00446 до 0.03889, увеличивается при нагреве внешнего или внутреннего слоев. Локальные и глобальные экстремумы распределения излучательной способности ЭВТИ во времени соответствуют локальным и глобальным экстремумам противоположного характера теплового сопротивления ЭВТИ – минимумы соответствуют максимумам и наоборот.

4. Термическое сопротивление ЭВТИ изменяется от 6.2 до 55.6 м² К/Вт и зависит от величины и направления теплового потока.

5. Термическое сопротивление сотовой панели изменяется от 0.41 до 9.4 м² К/Вт и имеет точку разрыва на границе режимов охлаждения криоэкранами внешней поверхности ЭВТИ и нагрева внутренней поверхности сотовой панели.

6. Наложение на распределения экспериментальных температур дополнительного малого возмущения, примерно равного экспериментальному, не нарушает характер идентифицированных функций. Глобальные экстремумы функций сохраняются приблизительно в том же временнóм диапазоне, за исключением сопротивления сотовой панели, глобальный максимум которого переходит из режима охлаждения “Д4” в аналогичный режим “Д2”.

Работа выполнена в рамках госзадания № 9.9074.2017/БЧ Минобрнауки России.