ISSN 0868-5886 НАУЧНОЕ ПРИБОРОСТРОЕНИЕ, 2022, том 32, № 4, c. 107-123
МАТЕМАТИЧЕСКИЕ МЕТОДЫ И МОДЕЛИРОВАНИЕ -В ПРИБОРОСТРОЕНИИ
УДК 536.2.083 © Н. О. Борщев, 2022
ПАРАМЕТРИЧЕСКАЯ ИДЕНТИФИКАЦИЯ КОЭФФИЦИЕНТА ТЕПЛОПРОВОДНОСТИ ЗЕРКАЛЬНЫХ МАТЕРИАЛОВ, ФУНКЦИОНИРУЮЩИХ НА ОРБИТАЛЬНОМ УЧАСТКЕ ПОЛЕТА
В данной работе представлен метод определения коэффициента теплопроводности зеркальных материалов космических аппаратов, функционирующих на орбитальном участке полета. Данная задача решается как задача поиска глобального экстремума параметризированного коэффициента теплопроводности в ходе минимизации среднеквадратичного функционала невязки между теоретическим и экспериментальным полями температур в местах установки датчиков температур. Обратные задачи считаются некорректными из-за "за-шумленных" входных данных, для преодоления этого необходимо применить регуляризацию. В данной работе используется метод итерационной регуляризации, где в качестве регуляризируемого параметра выступает номер итерации. В качестве алгоритма оптимизации выбран метод сопряженных градиентов как наиболее точный метод сходимости первого порядка.
Кл. сл.: обратная задача теплопроводности, метод итерационной регуляризации, среднеквадратичная ошибка, температурное поле, космический аппарат
ВВЕДЕНИЕ
При проектировании теплового облика изделий или оценке их термонапряженного состояния необходимо иметь априорную информацию об их исходных данных, таких как начально-граничные условия и функции изменения теплофизических параметров от факторов эксплуатационного применения. В данной работе рассматривается метод идентификации эффективного коэффициента теплопроводности материала как функции от температуры на примере элемента зеркальной системы космического аппарата (КА), который имеет сложную физико-химическую структуру и функционирует на орбитальном участке полета. При моделировании теплового режима данного объекта необходимо знать его эффективные теплофизи-ческие свойства, поскольку процесс теплообмена сопровождается как лучисто-кондуктивным теп-лопереносом на поверхности материала, так и переносом лучистой тепловой энергии внутри него.
ПОСТАНОВКА ЗАДАЧИ ТЕПЛООБМЕНА
Первым шагом в восстановлении коэффициента теплопроводности конструкции является составление "прямой" задачи теплообмена, моделирующей условия эксперимента.
Применительно к данной постановке задачи можно использовать упрощенную одномерную
схему теплопроводности при допущении о малом перераспределении теплового потока по поверхности конструкции при ее одномерном нагреве:
Сэф (T (x,r))
- А"
дх
8T (х,т)
дт
Кф (T(х,т))
дТ (х, т) дх
, дЧизл (х,т) . Т - Т .
~г _ 1 -*- -*- Г, 1
дх
дЧизл (х,т)
- е(Т(х,т))Лэф (Т(х,т))
дТ (1х ,т)
дх 4 4 " '' эф ^ '' дх
- Чизл(х,т)$(х,т)п(х).
Здесь:
Т — температура, К;
х — линейная координата, м;
т — время, с;
1х — длина образца;
Сэф(Т — эффективная объемная теплоемкость
Дж.
материала,
К
Кэф (Т) — эффективная теплопроводность
ма-
териала,
Вт . м • К
аизл (х,т) — удельный тепловой поток, погло-
Вт
щаемый оптическим зеркальным слоем, —-;
м
е(Т) — отражательно-излучательная характеристика материала;
п(х) — коэффициент преломления среды.
Выражение для толщины оптического слоя 3 имеет вид:
е(Т (х,т) xdx
^=![е ¡Л
Граничные условия имеют следующий вид: дТ (0,т)
X
дх ;
дТ(1х^=А (Та т)ч а
+£
(Т(I ,т))Го + а —аТ(I ,г)41,
\ V х' переизл 1 атм V х' ' I'
где
As (Т) — поглощательная способность материала в солнечном спектре;
— солнечный удельный тепловой поток,
Вт
м
датм — атмосферный удельный тепловой по-Вт
ток,
м
где:
¡их — относительный мидель поверхности, на которую падает солнечный поток, по направлению на Солнце;
а — среднее альбедо планеты;
50 — солнечная постоянная в окрестности планеты;
Сь С2 — константы, определяющие собственное излучение планеты, при этом для планет с плотной атмосферой (в частности, для Земли)
С _ (1— . с _ о.
4 ; 0;
ф1 — угловой коэффициент между поверхностью и планетой, определяющий долю собственного излучения планеты, попадающую на поверхность;
ф2 — комбинированный угловой коэффициент, зависящий от взаимного положения поверхности, Солнца и планеты, определяет долю отраженной от планеты солнечной энергии, попадающую на рассматриваемую поверхность.
а — постоянная Стефана - Больцмана, Вт4 ;
м К
Чжреизл — падающий переизлученный тепловой
Вт
поток между элементами конструкции, —-.
м
ОПРЕДЕЛЕНИЕ ПАДАЮЩЕГО УДЕЛЬНОГО ИНТЕГРАЛЬНОГО ТЕПЛОВОГО ПОТОКА
Плотности падающих тепловых потоков определяются по следующим формулам [1]: прямой солнечный поток
ас _ мА;
отраженный от планеты солнечный поток
аТ1 _
собственное излучение планеты _ С> + С2>
р
4 . Ч п
Н + Я
/ S e
Г /
В'
Г1-' Я ? 1
Рис. 1. К расчету угловых коэффициентов > и >
+
Векторы и углы, фигурирующие в формулах, приведены на рис. 1.
Относительный мидель поверхности определяется следующим выражением [1]:
(n, S) + |(n, S)|
2
где _п — единичный вектор нормали к поверхности; S — единичный вектор, направленный на Солнце.
Направляющие косинусы векторов п и S определяются в связанной с аппаратом правой системе координат.
Углы ц, у. и д. представляют собой:
ц — угол между нормалью к площадке F, для которой рассчитываются падающие лучистые потоки, и направлением на центр планеты;
у. — угол между направлением на площадку F из центра Земли и вектором S направления на Солнце;
д.. — угол между вертикальной плоскостью, содержащей нормаль к поверхности, и вертикальной плоскостью, содержащей вектор направления на Солнце.
Величина в0 (угол полуобзора планеты) определяется следующим образом:
в0 = arcsm I-I,
0 + Н )
где Н — высота полета.
Формула для определения углового коэффициента ф\ имеет вид:
cos^sin2 в0 при0 <у <--в0;
Pi =
cos у sin2 в0
ж
ж
— + arcsin(ctg у ctg в0
+
1
+ — п
Vsin2 во
cos у
sin у
— COSв0^Js\ñ2в'0 ж
cos у
точное решение в интервале 0 <ц < ж / 2 -в0 (поверхность F в этом случае не пересекает планету) и аппроксимируют точное решение с погрешностью менее \% при ж / 2 -в0 <ц < ж / 2 + А0, при этом функции выражаются следующим образом [\]:
/*(0о,ц)Я\(Ро,¥) = k АоМАо,ц);
sin2(eo)
/2(во) = \
л ■ 2п -> • зл cos в0 1 -sine0 1 + sin2 в0 + 2sinJe0 +—:—-ln 0
V
2sine 1 + sine
0 J
/3* (в0,у) =
Í /ъ (в0 )
/3 (в0 )-
^ ж
в0 +--V
2в
при 0 < у < п/2-в0,
при п/2-в0 <у < п / 2 + в0;
(в ) = cos2 в0(Ъ + sin2 в0) 1n 1 - sin в0
16sine
1 + sin в
ж „ ж п пРи--в0 < — + в0.
Угловой коэффициент ф2 вычисляется следующим образом:
P2 = /* (в>, у) cos(ys) + f* sin у sin(ys) cos S,,
где функции /2*(в0,у) и /3*(в0,у) определяют
(1 - sin0o)(3 + 3sinsin0o + 2sin2 в0) 8 '
Расчет углов проводится на основе траектор-ных параметров (перицентра, апоцентра орбиты, наклонения, аргумента перицентра, даты старта) и ориентации КА.
Угол ys находится из скалярного произведения векторов S и n0
cos ys = Sno;
n0 — единичный вектор нормали к плоскости орбиты КА, дополняет радиус-вектор КА и вектор скорости КА до правой системы координат.
Для ys после тригонометрических преобразований получаем следующее выражение:
sin ys = COs(i)sÍn(S) - sm(i)sÍn(y)cOs(S),
где i — наклонение орбиты (угол между плоскостью небесного экватора и плоскостью орбиты); y= v- Q (v — прямое восхождение Солнца: угловое расстояние по небесному экватору от точки весеннего равноденствия до меридиана Солнца, Q — долгота восходящего узла, отсчитываемая от точки весеннего равноденствия).
Величина 8 определяется из соотношения
8 = arcsm[sin(235)smv.
Атмосферный удельный поток равен сумме теплового потока за счет рекомбинации атомов и молекулярного потока:
Я-атм а + Ямол .
Удельный тепловой поток за счет рекомбинации атомов вычисляется по формуле [1]:
F
Я _ дЕ,Ни^, F
где N — число атомов газа в единице объема; С — эффективность рекомбинации атомов газа (С < 1); Е1 — энергия рекомбинации на один атом газа; и — скорость движения КА по орбите; — площадь проекции поверхности F на плоскость, перпендикулярную направлению скорости, Е — площадь поверхности F.
Удельный молекулярный тепловой поток равен:
Ямол _ 1 «Р(Н,
2
кг
Чп.
_ Чпад (1 - Si )фц ,
где япад — падающий удельный тепловой поток, Угловой элементарный коэффициент излуче-
м
ния фу может быть найден из выражения: cos(® ) cos(®y)
rn 2
зависимость искомых теплофизических характеристик от температуры. В данной работе используются линейно-непрерывные базисные функции, имеющие следующий вид [2-4]:
Nm (T) _\
0
T - T
Т - T
ПРИ Т < ^
при Tm-1 < Т < Tm ПРИ T > Tm,
m е 1, M;
параметрическое значение коэффициента теп-
лопроводности зеркального материала,
Вт м • К'
где а = 0.9^1.0 — коэффициент аккомодации; р —
плотность атмосферы на данной высоте, —3.
м
Выражение для переизлученного теплового потока имеет вид:
Определим количество M временных блоков, в каждом из которых Km (m е 1,M) шагов т по времени и на каждом из которых коэффициент Ар в линейных комбинациях постоянен.
Это количество определим из верхней оценки функциональной невязки:
K 82 < 8 ,
m f sum'
где Km — число узлов с замерами температуры по времени; Ssum — суммарная погрешность, Sf — ошибка замера температуры.
Отсюда получаем количество Km временных шагов т в каждом m-м блоке
K _ 8sum
8
индексы /, у — номера излучающих поверхностей элементов; ¡у — расстояние между центрами поверхностей, м; ю, у — углы между нормалью к излучающей поверхности и направлением на излучающую поверхность.
АЛГОРИТМ ИДЕНТИФИКАЦИИ ЭФФЕКТИВНОГО КОЭФФИЦИЕНТА ТЕПЛОПРОВОДНОСТИ
Представим теперь эффективную теплопроводность через коэффициенты в параметризирован-ном виде, умноженные на соответствующие базисные функции, учитывающие зависимость от температуры:
м
Хэф (Т) «X Хрнт (Т),
т_1
где Нт(Т) — базисные функции, описывающие
Если теперь весь временной промежуток [0, тшах] разделить на число Кт временных слоев в каждом блоке, получим количество М конечных элементов
т
м _-
K
а длина ДТт каждого конечного элемента т е 1,М равна:
Т — Т . -
ДТ __тах-пш, т е 1,М.
т м
Максимальное значение температуры может быть оценено из следующей постановки задачи [5-8]:
С* (Т (т) )дТТТ) _ АЛах(т)Яя +
+ ^Ттах (тТ[Чпереизл + Чатм — аТшт (Т) ];
Ттах (0) _ Т0.
0
А
p
Рис. 2. Схема теплового нагружения объекта испытаний.
Звездочками схематично изображены места установки термопар
А минимальное значение температуры — из выражения:
Сэф (Т (г) )дТд^ = = кф (Т)Ттт(г) [Ттах(г) - Ттп(г)] + дизл (1х ,г);
Т^п(0) = То.
Решив эту систему, было определено оптимальное количество временных блоков М, равное 3. Везде далее суммирование ведется по ним.
Схематичный вид образца приведен на рис. 2.
Рассмотрим восстановление искомых тепло-физических характеристик прибора на основе среднеквадратичного функционала невязки между теоретическим и экспериментальным полями температур [9— \\]:
\ Гт„
5(к) = 2 I Е[Т(к) -Т(х,г)]2dг.
В работе используется метод безусловной минимизации функционала 5 (к) с помощью градиентного метода сопряженных направлений как наиболее точного метода первого порядка, позволяющего достичь требуемой сходимости за минимальное число итераций.
Последовательный алгоритм метода сопряженных градиентов можно представить в следующем виде [12—17]:
где
х ;+1 = х np +ах
АХ ;+1 =-рк pn
Направление спуска определяется из соотношений:
pn = grad(Xn) + Pnpn-1; P0 = 0, p0 = grad S (Х0);
I |2
|grad S (Хn )|
Pn i2 •
|grad S (Х n-1)|
Погрешность входных данных Ssum, вычисленная в той же метрике, что и целевой функционал:
S = S + S, + S + S ,
sum a f окр х '
где:
Sf — погрешность входных температур, определяемая следующим выражением:
'max M
Sf = J T^Sl(T)dr;
0 L=\
докр — погрешность округления полученных значений;
да — случайная погрешность измерений; дх — погрешность постановки задачи, определяемая как
S =-
Г (Т(х -))дт(х,т) д Сэф(т М) дт ~дх Wtv \\дТ (х,т) ^ (т (х,т)) дх J дЧизл (х,7) дх
д дх 1 m \\дТ (х,т) (т() дх J
Как видно из выбранного метода оптимизации,
первостепенной задачей является отыскание гра- д5(кр) = Т^ гТ к ) Т(__)1дТ(к
диента среднеквадратичной интегральной ошибки: дк ~ ^ ( р (х,Т)] дк Г.
илр 0 иЛр
Для нахождения производной температуры по параметризированному значению теплопроводности достаточно продифференцировать основную постановку задачи теплообмена:
дСэф (Т(х,т)) дТ(Т(х,т))
дТ
+ Сэф (Т(х,т))
дт
дв( х,т)
в( х,т) +
дт
д_ дх
дХэф (Т(х,т)) дТ(х,т)
дТ
+Хэф (Т (х,т))
дх
дв( х,т)
дх
в( х, т) +
+
д2Чизл (х,т)
дхдХ„
Производная от толщины оптического слоя:
д5(х,Т) _ 1 \ в(т)х дв(х,у) _
дТ ¡х 0 дТ
-xdx.
Граничные условия примут следующий вид: дв(0,т)
дх
_ 0;
дХэф(Т) дТ(¡х,т) _ дАя (Т(¡х,т)) дТ(¡х,т)
дТ
дх
дТ
дх
Яя +
+
де<Тт))дТ(¡x,т) ([а + я ] — аТ(¡,т)4) —
дТ дх переизл ¿з* V х' / у
5 (Хр +ДХр) _ 5 (Хр) — «
5 (Хр + ДХр) _
д® (Хр) дХ
2 х г
дТ X ) Т(Х ) — а ——^
V р) с дх
— Т (х,т)]^т.
Откуда, согласно принципу глобального минимума, необходимо и достаточно приравнять полученное выражение нулю и выразить шаг спуска.
Получим:
а
У Т(Хр) — Т(х,т) ] —^-dт.
дТ (Хр) дХ
Таким образом, данный алгоритм включает в себя следующую последовательность действий:
1. Задание начальных приближений параметрических величин теплофизических характеристик.
2. Решение "прямой" задачи прогрева конструкции в одномерном приближении, моделируя условия наземной тепловакуумной отработки изделия.
—4£(Т (¡х ,т))аТ (¡х ,т)3в(¡x ,т).
Данная постановка задачи решается конечно-разностным методом по явной схеме, где конечно-разностная аппроксимация температур в узлах уже известна из решения "прямой" задачи прогрева.
Для нахождения шага спуска асп, исходя из метода итерационной регуляризации [9], запишем выражение целевого функционала на следующей итерации:
Рис. 3. Блок-схема алгоритма идентификации искомых теплофизических параметров
2
2
3. Получение экспериментального температурного поля изделия в местах установки датчиков температур.
4. Составление среднеквадратичной интегральной ошибки между теоретическим и экспериментальным температурными полями в местах установки датчиков температур.
5. Решение двух сопряженных задач по поиску компонент градиента целевого функционала невязки между теоретическим и экспериментальным температурными полями.
6. Вычисление шага спуска в методе сопряженных направлений на основе метода итерационной регуляризации.
7. Получение следующих итерированных приближений искомых параметрических величин.
8. Проверка критерия останова итерационного процесса. В случае его выполнения параметризи-рованные величины считаются искомыми, иначе — необходимо повторно выполнить пункты 1-7.
Реализация решения данного алгоритма проиллюстрирована блок-схемой на рис. 3.
ЧИСЛЕННЫЙ ВИРТУАЛЬНЫЙ ЭКСПЕРИМЕНТ И РЕЗУЛЬТАТЫ РАСЧЕТА
В численном эксперименте рассматривается образец в виде параллелепипеда, по толщине которого установлено 6 термопар (рис. 2). Все его поверхности теплоизолированы, кроме верхнего основания, на которое падает лучистый интегральный тепловой поток. Таким образом, реализуется одномерный прогрев по толщине материала, моделирующий заданную постановку задачи. В качестве источника теплового потока используется медный линейчатый нагреватель. Результаты определения падающего теплового потока и экспериментального температурного поля представлены на рис. 4 и 5.
1
2 \
ч
V. -
-с[с]2ст+с1рсгс|2] - Ам]\(Чп
20
40 60
Время, мин
80
100
Рис. 4. Удельные падающие тепловые потоки на нагреваемую поверхность образца. 1 — солнечный тепловой поток, 2— инфракрасный земной и атмосферный тепловые потоки
При итерационном уточнении параметризиро-ванной величины коэффициента теплопроводности материала будет также по итерациям восстанавливаться температурное поле, стремясь к своему экспериментальному аналогу. На рис. 6 приведены кривые температур в точках замера в разные моменты времени в зависимости от номера итерации.
Как видно из рисунка, для итерационной сходимости к своему итерационному постоянному значению необходимо 3 итерации, что говорит об эффективности предложенного метода.
Нагляднее всего процесс сходимости показан по минимизации среднеквадратичного отклонения теоретического температурного поля от экспериментального в местах замера температур. Данный процесс показан на рис. 7.
Рис. 6. Итерационное изменение температурного поля
Рис. 7. Итерационное изменение среднеквадратичной ошибки между теоретическим и экспериментальным температурными полями в местах установки датчиков температур.
1-3 — временные блоки
Итерационные изменения коэффициента теплопроводности материала представлены для каждого из временных блоков на рис. 8, а изменение коэффициента теплопроводности от температуры показано на рис. 9.
ВЫВОДЫ
1. Разработан метод параметрической идентификации коэффициента теплопроводности зеркальных материалов как функции от температуры методом итерационной регуляризации в приближении однонаправленного прогрева для КА на орбитальном участке полета.
2. Продемонстрированы результаты данного алгоритма на примере определения коэффициента теплопроводности элемента зеркальной поверхности космической обсерватории.
3. Результаты показали, что при данном уровне температур коэффициент теплопроводности будет
Вт
варьироваться в пределах 16.5-21 -.
м • К
4. Данный алгоритм может быть использован и для более широкого температурного диапазона для определения коэффициента теплопроводности материалов в изотропном приближении.
СПИСОК ЛИТЕРАТУРЫ
1. Залетаев В.М., Капинос Ю.В., Сургучев О.В. Расчет теплообмена космического аппарата. М.: Машиностроение, 1979. 208 с.
2. Крейн С.Г., Прозоровская О.И. Аналитические полугруппы и некорректные задачи для эволюционных уравнений // Доклады Академии наук СССР. 1960. Т. 133, № 2. С. 277-280.
URL: http://mi.mathnet.ru/dan23812
3. Басистое Ю.А., Яновский Ю.Г. Некорректные задачи в механике (реологии) вязкоупругих сред и их регуляризация // Механика композиционных материалов и конструкций. 2010. Т. 16, № 1. С. 117-143. URL: https://www.elibrary.ru/item.asp?id=15056455
4. Бакушинский А.Б., Кокурин М.Ю., Кокурин М.М. Прямые и обратные теоремы для итерационных методов решения нерегулярных операторных уравнений и разностных методов решения некорректных задач Ко-ши // Журнал вычислительной математики и математической физики. 2020. Т. 60, № 6. С. 939-962. DOI: 10.31857/S0044466920060022
5. Ефанов В.В., Мартынов М.Б., Карчаев Х.Ж. Летательные аппараты НПО им. С.А. Лавочкина (к 80-летию предприятия) // Вестник НПО им. С.А. Лавочкина. 2017. № 2. С. 5-16.
URL: https://www.elibrary.ru/item.asp?id=29237867
6. Блох А.Г., Журавлев Ю.А., Рыжков Л. Н. Теплообмен излучением. Справочник. М.: Энергоатомиздат, 1991. 432 с.
7. Тулин Д.В., Финченко В.С. Теоретико-экспериментальные методы проектирования систем обеспечения теплового режима космических аппаратов // Проектирование автоматических космических аппаратов для фундаментальных научных исследований. М.: МАИ-ПРИНТ. 2014. Т. 3. с. 1320-1437. URL: https://www.elibrary.ru/item.asp?id=25577466
8. Цаплин С.В., Болычев С.А., Романов А.Е. Теплообмен в космосе. Самара: изд-во "Самарский университет", 2013. 53 с. URL: http://tf.samsu.ru/pdf/heat_in_space.pdf
9. Алифанов О.М., Артюхин Е.А., Румянцев С.В. Экстремальные методы решения некорректных задач. М.: Наука, 1988. 288 с.
10. Алифанов О.М. Обратные задачи теплообмена. М.: Машиностроение, 1988. 280 с.
11. Формалев В.Ф. Теплоперенос в анизотропных твердых телах. М.: Физматлит, 2015. 238 с.
12. Васин В.В. Модифицированный метод наискорейшего спуска для нелинейных регулярных операторных
уравнений // Доклады Академии наук. 2015. Т. 462, № 3. С. 264. DOI: 10.7868/S0869565215150086
13. Голичев И.И. Модифицированный градиентный метод наискорейшего спуска решения линеаризованной задачи для нестационарных уравнений Навье-Стокса // Уфимский математический журнал. 2013. Т. 5, № 4. С. 60-76. URL:
https://www.elibrary.ru/item.asp?id=20930477
14. Формалев В.Ф., Ревизников Д.Л. Численные методы. М.: Физматлит, 2004. 400 с.
15. Формалев В.Ф. Анализ двумерных температурных полей в анизотропных телах с учетом подвижных границ и большой степени анизотропии // Теплофизика высоких температур. 1990. Т. 28, № 4. С. 715-721.
16. Формалев В.Ф. Идентификация двумерных тепловых потоков в анизотропных телах сложной формы // Инженерно-физический журнал. 1989. Т. 56, № 3. С. 382386.
17. Формалев В.Ф., Колесник С.А. Аналитическое решение второй начально-краевой задачи анизотропной теплопроводности // Математическое моделирование. 2003. Т. 15, № 6. С. 107-110.
Астрокосмический центр ФГБУ "Институт им. С.А. Лебедева ", Москва
Контакты: Борщев Никита Олегович, [email protected]
Материал поступил в редакцию 05.08.2022
ISSN0868-5886 NAUCHNOE PRIBOROSTROENIE, 2022, Vol. 32, No. 4, pp. 107-123
PARAMETRIC IDENTIFICATION OF THE THERMAL CONDUCTIVITY COEFFICIENT OF MIRROR MATERIALS OPERATING IN THE ORBITAL FLIGHT SECTION
N. O. Borshchev
Astrocosmic Center of the Federal State Institution of Science S.A. Lebedev Institute, Moscow, Russia
This paper presents a method for determining the thermal conductivity coefficient on mirror materials of spacecraft operating in the orbital flight section. This problem is solved as a problem of searching for the global extremum of the parametrized coefficient of thermal conductivity in the course of minimizing the root-mean-square functional of the discrepancy between the theoretical and experimental temperature fields in the places where temperature sensors are installed. Inverse problems are considered incorrect due to "noisy" input data. To overcome this, it is necessary to apply regularization. This study employs the iterative regularization method, with the iteration number serving as the regularized parameter. The method of conjugate gradients is chosen as the optimization algorithm, because it is the most accurate method of the first order of convergence.
Keywords: inverse heat conduction problem, iterative regularization method, root-mean-square error, temperature field, spacecraft
INTRODUCTION
When designing the thermal appearance of products or assessing their thermal stress states, it is necessary to have a priori information about their initial data, such as initial boundary conditions and dependence of changes in thermophysical parameters on operational application factors. In this paper, we consider a method for identifying the effective thermal conductivity of a material as a function of temperature using the example of a spacecraft (SC) mirror system component, which has a complex physical and chemical structure and operates in the orbital flight segment. When modeling the thermal mode of this object, it is necessary to know its effective thermophysical properties, since the heat transfer process is accompanied by conduction-radiation heat transfer both on the surface of the material and the transfer of radiant thermal energy inside it.
STATEMENT OF THE HEAT TRANSFER PROBLEM
The first step in defining the thermal conductivity of the structure is the determination of a "direct" heat transfer problem that simulates the conditions of the experiment.
With regard to this formulation of the problem, a simplified one-dimensional heat conduction scheme can be used under the assumption of a small redistribution of the heat flow over the surface of the structure during its one-dimensional heating:
C3é (T ( x,z))
ôT ( x,T)
ôz
ôx
k (T(x,z))
ôT ( x,z) ôx
, ^(x,z). T=T
ôx
ôq ( x,z)
±U3Jl V * '
ôx
= s(T(x,z))k3t (T(x,z))
ôT (lx ,z)
ôx
- qu3u( x,z)S(x,T )n( x).
Here:
T — temperature, K; x — linear coordinate, m; t — time, s; lx — the sample length;
C^T) is the effective volumetric heat capacity of J
the material, —, K
k (T) is the effective thermal conductivity of the W
material,
m • K
qU3m(x,t) is the specific heat flow absorbed by the W
optical mirror layer, —-;
m
s(T) is the reflective-radiative characteristic of
the material;
n(x) is the refractive index of the medium.
The expression for the thickness of the optical layer S has the form:
s(T ( x,r) xdx
5 = -Me
Li
The boundary conditions have the following form: 8t (0,t)
X
dx '
ST (lx ,t)
dx
= As (T(lx,0)qs
+s
(T(l ,t))l~q + q -oT(l ,r)41,
\ v x5 /1 nepeu3Ji ± amM v x ' ^ I '
where
^^ (r) is an absorption capacity of the material in the solar spectrum;
W
qs is the solar specific heat flow, —-;
m
W
qamM is the atmospheric specific heat flow, —-;
m
a is the Stefan - Boltzmann constant
W
' m2K4 '
qnepeu3J is incident reradiated heat flow between W
structural elements, —-;
m2
DETERMINATION OF INCIDENT SPECIFIC INTEGRAL HEAT FLOW
The densities of the incident heat flows are determined by the following formulas [1]: direct sunlight
qc _ mA;
solar flow reflected off the planet
qZL _
the planet's own radiation
qea6 _ ClVl + C2^2,
where:
is the relative midsection of the surface, on which the solar flow falls, in the direction of the Sun; a is the average albedo of the planet; So is the solar constant in the vicinity of the planet;
C1, C2 are the constants that determine the planet's own radiation, while for planets with dense atmospheres (in particular, for the Earth)
C _ (1 ~ a)So . c _ o-
C1 4 ; C2 o;
y1 is the angular coefficient between the surface and the planet, which determines the proportion of the planet's own radiation that hits the surface;
is the combined angular coefficient, which depends on the relative position of the surface, the Sun, and the planet and determines the fraction of solar energy that is reflected from the planet and hits the surface under consideration.
The vectors and angles appearing in the formulas are shown in Fig. 1.
The relative midsection of the surface is determined by the following expression [1]:
(n, S) + |(n, S)|
2
where n is the unit vector of the normal to the surface; S is the unit vector directed to the Sun.
The directional cosines of the vectors n and S are determined in the right coordinate system associated with the SC.
The angles y, ys and Ss are:
y is the angle formed by the normal to the site F, for which the incident radiant flows are calculated, and the direction to the center of the planet;
ys is the angle formed by the direction to the site F from the center of the Earth and the vector S directed to the Sun;
âs is the angle formed by the vertical plane containing the surface normal and the vertical plane containing the vector directed to the Sun;
The value of 60 (the planet's half-view angle) is determined as follows:
A • ( R )
0O = arcsin I-I,
0 U + H )
where H is the flight altitude.
The formula for determining the angular coefficient ^ is:
+
Pi =
2 л
cos у sin в0 при 0 <у <--в0;
cos y sin2 60
л
— + arcsin(ctg у ctg в0
i
+ — л
-y/sin2 во
cos y
Sin у
—cosв0^Jsin2 в0
л
cos y
лл
пРи — - в0 <у < — + в0.
The angular coefficient is calculated as follows: <P2 = /2* (0q , W) cos<7. ) + sin W sin<7. )cos Ss,
where the functions /*(0Q,y) and /*(0Q,y) determine the exact solution in the interval Q< y < n / 2- 0Q (the surface F does not intersect the planet in this case) and approximate the exact solution with an error of less than 1% for n / 2 -dQ <y < n / 2 + dQ. The functions are expressed as follows [1]:
/2*(0O,w) = pI(0Q,w) = k (0Q)p:(0Q,w);
sm2^)
/2(в0) = i
i + sin2 в + 2sin3вn +
cos4 вп, i - sin в ln
2sinв() i + sinв() j
/3 (в0,у) = Í /3 (в0 )
при 0 < у < л/2-в0
я л
/3в)-^- при л/2-в0 <у < Л/2 + в0;
/з(в0) =
2в
cos2 в0 (3 + sin2 в0) i - sin в
^тв0
ln-
i + sin в0
After trigonometric transformations, we obtain the following expression for ys:
sin ys = cos(i)sin(8) - sin(i)sin(y)cos(8),
where i is the inclination of the orbit (the angle between the plane of the celestial equator and the plane of the orbit); y= v- Q (v is the right ascension of the Sun: the angular distance along the celestial equator from the vernal equinox to the meridian of the Sun, Q is the longitude of the ascending point, counted from the vernal equinox).
The value of 8 is determined from the relation
8 = arcsin[sin(235)sinv.
The atmospheric specific flow is equal to the sum of the heat flow due to the recombination of atoms and molecular flow:
QamM q + Quoji .
The specific heat flow due to the recombination of atoms is calculated by the formula [1]:
F
q = gE, Nu-F,
where N is the number of gas atoms per volume unit; ^is the efficiency of recombination of gas atoms (£< < 1); Ei is the recombination energy per gas atom; u is the velocity of the SC in orbit; FN is the area of surface F projection onto the plane perpendicular to the velocity direction, F is the area of the surface F.
The specific molecular heat flow is equal to:
qMoj =1 aP(H )v3~F, 2 fn
where a = 0.9^1.0 is the accommodation coefficient; p is the density of the atmosphere at a given height, kg m3.
The expression for the reradiated heat flow has the form:
^переизл = Чпад (i - Sj
(1 - sin0o)(3 + 3sinsin0o + 2sin2 90) 8 .
Angles are calculated based on trajectory parameters (periapsis, orbit apoapsis, inclination, periapsis argument, launch date) and the orientation of the SC.
The angle ys s is found from the scalar product of the vectors S and n0
cos ys = Sno;
W
where qnad is the incident specific heat flow, —-. The
m
expression for angular elementary emissivity yij is:
Pu =-
cos(®i ) cos(®.. )
"Л2 :
indices i, j are the numbers of radiating surfaces of elements; lij is the distance between the centers of the
• •. . n .. , . cr^ u-4. surfaces, м; rn, . are the angles formed by the normal
n0 is the unit vector of the normal to the SC orbit , ,, } г Л, i- : i- .
plane, which complements the SC radius vector and the SC velocity vector to the right coordinate system.
to the radiating surface and the direction to the radiating surface.
+
IDENTIFICATION ALGORITHM FOR THERMAL CONDUCTIVITY EFFECTIVE COEFFICIENT
Let us now represent the effective thermal conductivity in terms of the coefficients in a parametrized form and multiplied by the corresponding basis functions, taking into account the dependence on temperature:
^ (T ) (T ),
where Nm(T) are the basis functions that describe the temperature dependence of the desired thermophysical characteristics. In this work, linearly continuous basis functions are used that have the following form [2-4]:
Nm (T ) = \
0
T - T
T_ - T
If T < Tm-1,
If Tm-1 < T < Tm
if T > T
m e 1, M ;
Xp is the parametric value of the thermal conductivity
of the mirror material,
W
m • K
Let us determine the number M of time blocks, in
6Sf'
K.,
and the length ATm of each finite element m e 1,M is:
T - T -
AT = -m e 1,M.
m M
The maximum temperature value can be estimated from the following problem statement [5-8]:
C* (T (t) = AsTmx(z)qs
+
+ STmax(T)[qnePeu3n + qamM -^Tmin(T)4];
Tmax(0) = To.
And the minimum temperature is taken from the expression:
c* (T (t) )
^Tmin(T) ÔT
each of which Km (m e 1,M) steps, lasting for t and at each the coefficient XP is constant in linear combinations.
This number is determined from the upper estimate of the functional residual error
K S2 < S ,
m f sum 5
where Km is the number of sites with temperature measurements; Ssum is the total error, Sf is the temperature measurement error.
From here, we obtain the number Km of time steps t in each m-th block
K _ S sum
_ (T)Tmin T) [Tmax T) " Tmrn(T)] + ^ (lx ,t);
Tnin (0) _ T0.
Having solved this system, the optimal number of time blocks M was determined, equal to 3. Everywhere below, the summation is carried out using it. A schematic view of the sample is shown in Fig. 2.
Fig. 2. Scheme of thermal loading of the test object. The asterisks show the places where thermocouples are installed
Let us consider the determination of the desired thermophysical characteristics of the device based on the root-mean-square residual functional between the theoretical and experimental temperature fields [911]:
1 Tux
s(*P) _ 2 i E[t(*p) -T(x,f)]2dr.
2 0
The study uses the unconditional minimization of the functional S(Xp) using the gradient method of
conjugate directions as the most accurate first-order method that allows achieving the required convergence in the minimum number of iterations.
The sequential algorithm of the conjugate gradient method can be represented as follows [12-17]:
If we now divide the entire time interval [0, Tmax] by the number Km of time steps in each block, we get the number M of finite elements
T
M _-
where
Xn+1 = Xn + AXn+ p p p
ax ;+1 =-ßk pn.
The direction of descent is determined by the ratios:
p" = grad(X") + PnP"- ;
0
ß = 0, p0 = grad S (X0);
ßn =
I grad S (Xn )| Igrad S (X n"1)|:
The input data error Ssum, calculated using the same metric as the objective functional, is:
S = S + Sf + S + 8
sum a f окр x
where:
Sf is the input temperature error, defined by the following expression:
max M
Sf = i £8,(r)dr;
0 L=1
SOKp is the rounding error of the obtained values;
Sa is a random measurement error;
Sx is the problem statement error, defined as
8 =-
Г (T(г ,))ÔT(x,,) Ô C* (T ( x,,) ) Ô, "ôx Wtv \\ÔT ( x,) [Я (T(x,T)) Ôx \ дЧгзл (x,,) Ôx
ô Ôx Wtv \\ÔT ( x,) [Я (T(x,T)) Ôx \
As can be seen from the chosen optimization method, the primary task is to find the gradient of the mean square integral error:
ascAp)
дЯ„
j£[T (Я, ) - T ( XJ)fTk- d,
dk„
дСЭф (T ( x,,) ) ÔT (T ( x,,) )
ÔT
+ Сэф (T(x,,))
Ô, дв( x,,)
в( x,,) +
Ô,
Ô ГдЯэф (T(x,,)) ÔT(x,)
Ôx ÔT Ôx
, a ITÎ \\Ôe(x,)
+Яэф (T(x,,))
Ôx
в( x,) +
Ô"Чизл (x,,)
+
ÔxÔA
The derivative of the thickness of the optical layer:
ÔS( x,T ) = _1jrge(r ) x
Ôs( x,y)
ÔT lx 0 ÔT
xdx.
The boundary conditions take the following form: ô0(O,z)
дЯэф(T) ÔT(lx,)_ ÔAs (T(lx,)) ÔT(lx,,)
ÔT
Ôx
ÔT
q* +
+
ÔS(T(l■ ^)M,([q + q ]-aT(l ,,)4)
To find the derivative temperature according to the parameterized value of thermal conductivity, it is sufficient to differentiate the main formulation of the heat transfer problem:
Ôx
= 0;
-4e(T (lx ,r) )aT (lx ,r)'e(lx ,r).
This problem is solved by the finite-difference method according to an explicit scheme, where the finite-difference approximation of temperatures at the points is already known from the solution of the "direct" heating problem.
To find the descent step acn, based on the iterative regularization method [9], we write the expression for the objective functional for the next iteration:
S (Я, +ДА ) = S (Я, )-aß
S (Я, +ДЯ, ) =
2
£ Г
t (Я, ) -ае
ÔT (Я, )
дЯ„
ÔS (Я, )
дЯ„
- T(x,f)]2dr.
Based on this, according to the global minimum principle, it is necessary and sufficient to equate the resulting expression to zero and express the descent step.
We obtain:
=£J
t (Я, ) - T ( x,)
ÔT (Я, ) дЯ
d,
2
2
1
Thus, this algorithm includes the following sequence of actions:
1. Setting the initial approximations of the values of thermophysical characteristics.
2. Solution of the "direct" problem of heating the structure in a one-dimensional approximation, simulating the conditions of ground thermal vacuum testing of the product.
3. Obtaining an experimental temperature field of the product at the installation sites of temperature sensors.
4. Determination of the root-mean-square integral error between the theoretical and experimental temperature fields at the installation sites of temperature sensors.
5. Solving two conjugate problems of finding the components of the gradient of the target functional of the discrepancy between the theoretical and experimental temperature fields.
6. Calculation of the descent step in the method of conjugate directions based on the method of iterative regularization.
7. Obtaining the following iterated approximations of the required parametric values.
8. Checking the criterion for stopping the iterative process. If it is executed, the parameterized values are considered to be the required ones, otherwise, steps 17 must be repeated.
The implementation of the solution of this algorithm is illustrated by the block diagram in Fig. 3.
Fig. 3. Block diagram of the algorithm for identifying the required thermophysical parameters
NUMERICAL VIRTUAL EXPERIMENT AND CALCULATION RESULTS
In the numerical experiment, a sample in the form of a parallelepiped is considered, inside which 6 thermocouples are installed (Fig. 2). All its surfaces are thermally insulated, except for the upper base, on which the radiant integral heat flow falls. Thus, one-dimensional heating throughout the entire thickness of the material is performed, simulating the given problem. A copper bar heater is used as a source of heat flow. The results of determining the incident heat flow and the experimental temperature field are shown in Figs. 4 and 5.
Fig. 4. Specific incident heat flows on surface of the sample to be heated.
1 — solar heat flux, 2 — infrared terrestrial and atmospheric heat flows
Fig. 5. Temperature field in the places where temperature sensors are installed. 1-6 — thermocouple numbers
The temperature field, tending to its experimental counterpart, is also restored by iterative refinement of the parametrized value of the thermal conductivity of the material. Fig. 6 shows the dependence of temperatures on the iteration number at the temperature measurement point over time.
Fig. 6. Iterative change of the temperature field
As can be seen from the figure, iterative convergence to its iterative constant value requires 3 iterations, which indicates the effectiveness of the proposed method.
Most clearly, the process of convergence is shown by minimizing the root-mean-square deviation of the theoretical temperature field from the experimental one at the points of temperature measurement. This process is shown in Fig. 7.
Fig. 7. Iterative change in the root-mean-square error between the theoretical and experimental temperature fields at the installation sites of temperature sensors. 1-3 — time blocks
Iterative changes in the thermal conductivity coefficient of the material are presented for each of the time blocks in Fig. 8, and the change in the thermal conductivity coefficient according to the temperature — in Fig. 9.
Fig. 8. Iterative dependence of the thermal conductivity coefficient of the material on temperature. 1-3 — time blocks
Fig. 9. The dependence of the coefficient of thermal conductivity of the material on temperature
CONCLUSIONS
1. A method has been developed for the parametric identification of the thermal conductivity coefficient of mirror materials as a function of temperature by the iterative regularization method in the approximation
of unidirectional heating for a SC in the orbital flight segment.
2. The results of this algorithm use are demonstrated on the example of determining the thermal conductivity of a component of the mirror surface of a space observatory.
3. The results showed that at a given temperature level, the thermal conductivity coefficient will vary
W
within 16.5-21 ——.
m • K
4. This algorithm can be used for a wider temperature range to determine the thermal conductivity of materials in the isotropic approximation.
REFERENCES
1. Zaletaev V.M., Kapinos Yu.V., Surguchev O.V. Raschet teploobmena kosmicheskogo apparata [Spacecraft Heat Exchange Calculation]. Moscow, Mashinostroenie Publ., 1979. 208 p. (In Russ.).
2. Krein S.G., Prozorovskaya O.I. [Analytic semigroups and incorrect problems for evolutionary equations]. Doklady Akademii Nauk SSSR [Reports of the USSR Academy of Sciences], 1960, vol. 133, no. 2, pp. 277-280.
URL: http://mi.mathnet.ru/dan23812 (In Russ.).
3. Basistov Yu.A., Yanovsky Yu.G. [Ill-posed problems of mechanics (rheology) of viscoelastic media and theirs re-gularization]. Mekhanika kompozitsionnykh materialov i konstruktsii [Mechanics of composite materials and structures], 2010, vol. 16, no. 1, pp. 117-143. URL: https://www.elibrary.ru/item.asp?id=15056455 (In Russ.).
4. Bakushinskii A.B., Kokurin M.Y., Kokurin M.M. [Direct and converse theorems for iterative methods of solving irregular operator equations and finite difference methods for solving ill-posed Cauchy problems]. Zhurnal vychisli-tel'noi matematiki i matematicheskoi fiziki [Computational Mathematics and Mathematical Physics], 2020, vol. 60, no. 6, pp. 939-962. DOI: 10.31857/S0044466920060022 (In Russ.).
5. Efanov V.V., Martynov M.B., Karchaev Kh.Zh. [Flightborne vehicles by Lavochkin association (to the eightieth anniversary of Lavochkin association)]. Vestnik NPO im. S.A. Lavochkina [Bulletin of the NPO named after S.A. Lavochkin], 2017, no. 2, pp. 5-16.
URL: https://www.elibrary.ru/item.asp?id=29237867 (In Russ.).
6. Blokh A.G., Zhuravlev Yu.A., Ryzhkov L.N. Teploobmen izlucheniem. Spravochnik [Heat exchange by radiation. Reference book]. Moscow, Ehnergoatomizdat Publ., 1991. 432 p. (In Russ.).
7. Tulin D.V., Finchenko V.S. [Theoretical and experimental methods for designing systems for ensuring the thermal regime of spacecraft]. Proektirovanie avtomaticheskikh kosmicheskikh apparatov dlya fundamental'nykh nauch-nykh issledovanii [Design of automatic spacecraft for basic scientific research]. Moscow, MAI-PRINT, 2014, vol. 3, pp. 1320-1437. URL:
https://www.elibrary.ru/item.asp?id=25577466 (In Russ.).
8. Tsaplin S.V., Bolychev S.A., Romanov A.E. Teploobmen v kosmose [Heat exchange in space]. Samara, Samarskii universitet, 2013. 53 p.
URL: http://tf.samsu.ru/pdf/heat_in_space.pdf (In Russ.).
9. Alifanov O.M., Artyukhin E.A., Rumyantsev S.V. Ehkstremal'nye metody resheniya nekorrektnykh zadach [Extreme methods of solving incorrect backdrops]. Moscow, Nauka Publ., 1988. 288 p. (In Russ.).
10. Alifanov O.M. Obratnye zadachi teploobmena [Reverse heat exchange problems]. Moscow, Mashinostroenie Publ., 1988. 280 p. (In Russ.).
11. Formalev V.F. Teploperenos v anizotropnykh tverdykh te-lakh [Heat transfer in anisotropic solids]. Moscow, Fizmatlit Publ., 2015. 238 p. (In Russ.).
12. Vasin V.V. [Modified quick descent method for nonlinear regular operator equations]. Doklady Akademii nauk [Reports of the Academy of Sciences], 2015, vol. 462, no. 3, pp. 264. DOI: 10.7868/S0869565215150086 (In Russ.).
13. Golichev I.I. [Modified gradient fastest descent method for solving linearized non-stationary Navier-Stokes equations] . Ufimskii matematicheskii zhurnal [Ufa Mathematical Journal], 2013, vol. 5, no. 4, pp. 60-76. URL: https://www.elibrary.ru/item.asp?id=20930477 (In Russ.).
14. Formalev V.F., Reviznikov D.L. Chislennye metody [Numerical methods]. Moscow, Fizmatlit Publ., 2004. 400 p. (In Russ.).
15. Formalev V.F. [Analysis of two-dimensional temperature fields in anisotropic bodies taking into account mobile boundaries and a large degree of anisotropy]. Teplofizika vysokikh temperatur [Thermophysics of high temperatures], 1990, vol. 28, no. 4, pp. 715-721. (In Russ.).
16. Formalev V.F. [Identification of two-dimensional heat flows in complex anisotropic forms]. Inzhenerno-fizicheskii zhurnal [Journal of Engineering Physics and Thermophysics], 1989, vol. 56, no. 3, pp. 382-386. (In Russ.).
17. Formalev V.F., Kolesnik S.A. [Analytical solution of the second initial-boundary problem of anisotropic thermal conductivity]. Matematicheskoe modelirovanie [Mathematical Models and Computer Simulations], 2003, vol. 15, no. 6, pp. 107-110. (In Russ.).
Contacts: Borshchev Nikita Olegovich, Article received ЬУ the editorial office on 05 08-2022