Математическое моделирование теплообмена в космических аппаратах

Введение

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

Технические возможности наземных комплексов зачастую не позволяют имитиро­вать тепловые режимы КА в полном объеме, в том числе в нештатных ситуациях. Мате­матическое моделирование при наличии до­стоверной тепловой математической модели КА позволяет прогнозировать поведение его систем как при наземных, так и при летных испытаниях. Использование математического моделирования при тепловой отработке КА позволяет:

• сократить объем ТВИ до автономных ТВИ составных частей (СЧ) КА;
• уменьшить количество режимов и про­должительность тепловакуумных испытаний;
• подтвердить тепловые режимы КА по результатам ТВИ СЧ КА и теплового расчета КА на основе математических моделей тепло­обмена.

Расчетно-экспериментальный метод эф­фективен только при его системности, призна­ком которой является наличие обратной связи между экспериментальными и расчетными данными, позволяющей уточнять тепловые математические модели непосредственно по результатам ТВИ. В качестве такой связи мо­гут выступать обратные задачи теплообмена (рис. 1).

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

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

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

Задачи идентификации уравнений тепло­вого баланса с распределенными параметрами ставятся при исследовании свойств перспек­тивных для космической техники теплоза­щитных материалов: волокнистых высоко­пористых [6-8], пористых реагирующих [9], неметаллических ультрапористых сетчатых [10-12] и анизотропных материалов [13, 14], слоистых сред [15], стеклопластиков [16], ар­мированных композитов [17], керамических материалов [18], тонкостенных композитных конструкций [19], разрушающихся теплоза­щитных материалов [20, 21], композицион­ных структур с микрокомпонентами, имею­щими различные теплопроводности [22, 23]. С помощью решения задач идентификации рассчитана тепловая проводимость радиаци- онно-кондуктивного теплообмена высокопо­ристых анизотропных полупрозрачных для теплового излучения материалов [24, 25], эф­фективная теплопроводность конструкцион­ных, теплозащитных и теплоизоляционных материалов [26]. Кроме того, такие задачи ста­вятся для изучения теплообмена в процессе литья и последующего остывания заготовки [27], параметров материала в процессе термо­диффузии [28], при квазистационарном оплав­лении и прогреве образца кварцевой стеклоке­рамики [29].

Проведена диагностика внешнего тепло­вого воздействия на элементы конструкции ав­томатической межпланетной станции Марс-96 на этапе ее выведения после сброса головно­го обтекателя [30]. Теплоперенос в элементах системы станции описывается трехмерным уравнением теплопроводности в различных системах координат на примере исследования внешнего теплового воздействия. Предложен расчетно-экспериментальный подход к диа­гностике нестационарного теплообмена с ис­пользованием масштабной модели летатель­ного аппарата, выполненной из сплава ниобия ВН-2 в виде треугольного крыла с затуплен­ными передними кромками, препарированной датчиками теплового потока с аппаратурой, которая построена на методах обратных задач теплопроводности [31].

С начала 90-х гг. выпускаются россий­ские долгоресурсные конкурентоспособные КА связи нового поколения с крупногабарит­ным негерметичным приборным отсеком в форме параллелепипеда, представляющим сборную на основе плоских прямоугольных трехслойных сотовых панелей блочно-модуль­ную конструкцию, несущую теплонагружен­ную бортовую аппаратуру (БА) и выполня­ющую одновременно силовую, тепловую и защитную от ионизирующих излучений функ­ции [32]. Для обеспечения теплового режима БА применяется пассивная система терморегу­лирования (ПСТР) на базе оптических покры­тий и нерегулируемых низкотемпературных тепловых труб (НТТ) разной конфигурации и типов профилей в сочетании с системами элек­троподогрева. К данному классу КА предъ­являются повышенные требования точности проектных оценок тепловых режимов БА и теплового состояния элементов конструкции в реальных условиях эксплуатации на геоста­ционарной орбите (ГСО), в том числе при воз­никновении аварийных ситуаций [33].

Представлены физические и динами­ческие тепловые математические модели в распределенных и сосредоточенных параме­трах, вычислительные алгоритмы, программ­ное обеспечение и некоторые результаты чис­ленных расчетов радиационно-кондуктивного теплообмена в несущих теплонагруженную бортовую аппаратуру блоках и модулях негер­метичного приборного отсека перспективных долгоресурсных КА связи с улучшенными по сравнению с традиционными герметичными отсеками массогабаритными и энергетически­ми характеристиками [32, 33].

Получены результаты расчетов по ма­тематическим моделям в сосредоточенных и распределенно-сосредоточенных параметрах при функционировании приборов бортовой аппаратуры и сети нерегулируемых тепловых труб в нетрадиционном Н-образном блоке мо­дуля служебных систем [34] и процессов радиационно-кондуктивного теплообмена мо­дуля полезной нагрузки [35] негерметичного приборного отсека перспективных долгоре­сурсных геостационарных КА негерметично­го блочно-модульного исполнения в условиях орбитальной эксплуатации.

Рассмотрены математические модели, учитывающие основные механизмы переноса тепла в приборных отсеках современных КА. Приведены результаты численного моделиро­вания теплофизических процессов в типичном модуле приборного отсека космического аппа­рата связи. На основе анализа многомерных температурных полей сделаны выводы о на­личии существенной неоднородности в рас­пределении температурного поля по приборам радиоэлектронной бортовой аппаратуры [36].

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

В настоящей работе продолжено форми­рование расчетно-экспериментального мето­да тепловой отработки КА с использованием обратных задач теплообмена. Объект иссле­дования вычислительных экспериментов - математическая модель с сосредоточенными параметрами негерметичного приборного от­сека КА.

Моделирование теплообмена в космических аппаратах

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

Как правило, уравнения взаимодействия включают функцию процесса и дифференци­альный оператор, характеризующий распре­деление потенциала внутри системы, и имеют первый (уравнение теплопроводности) или второй (волновое уравнение) порядок произ­водной по времени, а также второй порядок производной по координате.

Математические модели, основанные на дифференциальных уравнениях в частных производных, называются моделями с распре­деленными параметрами (МРП) (рис. 2, а). Если система приводится к конечному числу узлов, внутренние связи которых не рассма­триваются, для каждого узла записывается уравнение взаимодействия, в котором дифференциальный оператор по координатам равен нулю. Математические модели таких систем представляют собой системы обыкновенных дифференциальных уравнений и называют­ся моделями с сосредоточенными параметра­ми (МСП) (рис. 2, б). Взаимодействие узлов системы в МСП характеризуется коэффици­ентами взаимодействия, которые можно рас­считывать аналитически или эмпирически и уточнить экспериментально, а не свойствами непрерывной среды, как в МРП.

Для термодинамически открытой систе­мы в процессе первого порядка производной по времени справедливы следующие утвержде­ния [38]:

• коэффициенты взаимодействия связей внутренних элементов являются линейно за­висимыми;
• полный набор коэффициентов взаимо­действия не может быть определен только по граничным условиям и соответствующему им полю потенциала;
• при наличии приближения (l - 1) коэф­фициентов взаимодействия обратная задача идентификации коэффициентов последова­тельно решается в приближении l для (N- 1) связи каждого узла i системы;
• при отсутствии первого приближения полный набор внутренних связей восстанав­ливается по граничным условиям и полю по­тенциала только для (N - 1) связи одного узла.

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

Тепловые расчеты на основе математических моделей теплообмена в космических аппаратах

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

Целями проектных расчетов ТВИ СЧ КА являются разработка программы и теорети­ческое обоснование тепловых режимов ТВИ, по которым выполняется проверка тепловых математических моделеИ СЧ.

Цель проектных расчетов летных испы­тании СЧ КА заключается в теоретическом обосновании тепловых режимов летноИ экс­плуатации и выработке рекомендации по до­работке СЧ.

Цель поверочных тепловых расчетов ТВИ СЧ КА - проверка тепловых моделеи СЧ КА, идентифицированных по результатам проектных расчетов, при граничных услови­ях, соответствующих измеренным внешним и внутренним воздействиям в автономных ТВИ СЧ.

Цель поверочных тепловых расчетов лет­ных испытании СЧ КА - подтверждение те­пловых режимов СЧ, согласно предъявляемым требованиям для условии летнои эксплуатации на основе идентифицированнои тепловои математическои модели с точностью, установ- леннои в поверочных расчетах ТВИ СЧ.

Оценка точности тепловой модели с рас­пределенными параметрами осуществляется путем сопоставления значений температур в местах установки датчиков с показаниями датчиков, полученными в ТВИ. Точность теп­ловой модели выражается абсолютной  и относительной δΤ погрешностями темпера­туры, рассчитанными по контрольной выбор­ке результатов с установленными значениями доверительной вероятности, одна из которых равна надежности СТР СЧ. Абсолютная по­грешность температуры вычисляется в каждой контрольной точке n по температуре T_n,0 теп­ловой модели и показанию датчика Т_{n, 1} для каждого режима ТВИ [39]

В качестве истинного значения темпе­ратуры принято показание датчика. Как след­ствие, погрешности тепловой модели имеют смысл разностей теоретических и эксперимен­тальных температур и не включают ошибки измерения температуры в условиях ТВИ.

Абсолютная погрешность температуры СЧ определяется средним квадратичным для N точек контрольной выборки [39]

Относительная погрешность температу­ры СЧ равна частному среднего квадратично­го отклонения и допустимого диапазона тем­пературы

где T_min, T_max - границы диапазона.

Для оценки расхождения расчетных и экспериментальных данных в нестационарных режимах сформулировано предположение об аналогии погрешности расчета и функции от­клика электронного измерительного устрой­ства. Для процесса первого порядка производ­ной по времени функция отклика [39] примет следующий вид:

где у(τ) - функция отклика; τ - время (τ ≥ τ₀), с;

у - стационарное значение сигнала;

τ₀ - время подачи сигнала, с;

τ_0,63 - постоянная времени, с.

В момент времени (τ₀ + τ_0,63) значение функции составляет приблизительно 63 % ста­ционарного значения, а через 5τ_0,63 - 99 % (рис. 3).

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

Модель с сосредоточенными параметра­ми данной системы для процесса нагрева или охлаждения при наличии тепловой проводи­мости α через теплопередающее основание запишем в виде [39]:

где Τ - температура нагревателя; m - масса нагревателя, кг; с - удельная теплоемкость нагревателя, Дж/(кгК);

T_s - температура панели под основанием нагревателя, K.

Начальную температуру нагревателя обозначим T₀. Тогда

Решение полученного дифференциального уравнения следующее [39]:

Через бесконечный промежуток времени, соответствующий стационарному состоянию

где с_а - абсолютная теплоемкость гранично­го узла МСП, Дж/K;

α_α - абсолютная тепловая проводимость контакта граничного узла с одним из внутрен­них узлов, Вт/К.

Величина τ_0,63 является постоянной вре­мени рассматриваемой системы в процессе нагрева или охлаждения. Ее не сложно уста­новить экспериментально с помощью полу­чаемой экспоненциальной зависимости (см. рис. 3), поэтому частное абсолютной теплоемкости узла системы и тепловой проводимости контакта позволяет идентифицировать один из этих параметров.

Итерационная регуляризация задач идентификации

Для проверки модифицированного вариаци­онного метода итерационной регуляризации [40] проведен вычислительный эксперимент на примере идентификации тепловых прово­димостей математической модели составной части КА (рис. 4).

Идентификация выполнена на основе экспериментальных данных, полученных в тепловых вакуумных испытаниях составной части КА в режиме предельного нагрева от инфракрасных имитаторов плотностью тепло­вого потока до 900 Вт/м². В случае идентифи­кации проводимостей тепловых связей МСП прямым аналитическим методом определитель матрицы при векторе искомых проводимостей равен нулю. Следовательно, полный набор проводимостей тепловых связей неопределим по полям температуры и ее производной по времени, соответствующих заданным гранич­ным условиям. Однако если исключить первое уравнение и считать одну из проводимостей известной, при отсутствии первого прибли­жения система имеет единственное решение относительно неизвестных проводимостей тепловых связей первого узла с остальными узлами [38].

Для решения задачи с использованием МИР начальные приближения определены до­статочно произвольно. Дисперсия эксперимен­тальной температуры T_i ⁽⁰⁾ равна [39]

где M - количество точек временн0й сетки;

T′_ij - наиболее близкое к истинному значе­ние температуры, чаще всего среднее арифме­тическое ряда измерений в этой точке.

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

Следовательно, квадрат ошибки температур­ных измерений для всей системы равен [39]

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

Расчеты выполняются МИР [41] и моди­фицированным вариационным методом Ти­хонова [40] с определением параметра регу­ляризации по комбинированной методике для задачи идентификации с неустойчивым реше­нием (рис. 6, 7, табл. 1, 2). Единственность решения обеспечивается упрощением исход­ной системы до вида, в котором все искомые функции относятся к одному узлу и выража­ются через заданную в начальном приближе­нии функцию, которая и регуляризуется [38].

Устойчивость решений по Ляпунову про­веряется сопоставлением функций, определен­ных по возмущенному температурному полю, с идентифицированными функциями, которые считаются «истинными», а рассчитанное по ним температурное поле - «эксперименталь­ным». На «экспериментальное» поле наклады­вается нормальное возмущение (рис. 8).

Результаты применения метода итераци­онной регуляризации в случае нормального возмущения показаны на рис. 9, 10. В отли­чие от случая исходных экспериментальных температур идентифицированные функции определяются на итерации 40, а не на итера­ции 36, однако являются близкими по значе­ниям и характеру к идентифицированным при отсутствии возмущения. По этой причине их отличие от «истинных» функций, идентифици­рованных методом Тихонова, вызвано не нало­жением возмущений на «экспериментальное» температурное поле, а особенностями алгорит­мов самих методов.

Результаты применения вариационного метода Тихонова показаны на рис. 11, 12. Как и в случае исходных экспериментальных температур, идентифицированные функции опре­деляются на итерации 24 по условию замедле­ния сходимости температурного функционала и являются близкими по значениям и характе­ру к «истинным» функциям.

Применение математических моделей для тепловой отработки космических аппаратов

По идентифицированным параметрам по­строена математическая модель составной части КА, использованная при его тепловой отработке (рис. 13). Для разработки проект­ных и поверочных расчетов режимов тепло­вых вакуумных испытаний и летной эксплу­атации применялись идентифицированные МСП и МРП СЧ, обеспечившие погрешно­сти стационарных расчетов не более 9% от ширины допустимых диапазонов температу­ры для всех элементов СЧ за исключением экранно-вакуумной тепловой изоляции с относительной погрешностью 15% (табл. 3), (рис. 14).

Идентифицированные по двум группам датчиков функции времени излучательной спо­собности экранно-вакуумной теплоизоляции (ЭВТИ) относятся к диапазону 0,02...0,05 с доверительной вероятностью P = 1 - 3/162 ≈ 0,981. Среднее арифметическое диапазона равно 0,035. Отсюда следует, что излучательная способность ЭВТИ как системы из множе­ства слоев с коэффициентами излучения 0,04 и 0,05 проявляется как приведенный коэффи­циент излучения, получаемый в результате взаимодействия внутренних слоев и имеющий значение не выше максимального значения из- лучательных способностей этих слоев, а не как излучательная способность наружного слоя, равная 0,62.

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

Тепловые режимы СЧ КА подтверждены тепловым расчетом на основе математической модели теплообмена, идентифицированной по результатам ТВИ.

Выводы

1. Реализован метод тепловой отработки КА, основанный на тепловых расчетах и тепловых вакуумных испытаниях, с использованием обратных задач теплообмена для идентифи­кации математических моделей теплообмена.
2. Системный подход применен при на­земной тепловой отработке составной части космического аппарата в ОАО «Корпорация «Комета».
3. Применение идентифицированных параметров в тепловой математической моде­ли составной части КА позволило обеспечить погрешности стационарных расчетов не более 9 % от ширины допустимых диапазонов температуры для всех элементов за исключением экранно-вакуумной тепловой изоляции с отно­сительной погрешностью 15 %.