CC BY
78
Программирование методов расчета картины растекания токов по поверхности космических аппаратов
Борисов Н.И., Востриков А.В., Жадов А.Д., Колтунова В.Э.
Национальный исследовательский университет Высшая школа экономики borisovÇcù,itas.miem.edu.ru, avostrikov&hse.ru, exfaust(a)yandex.ru, harbin.lera(cp,yandex.ru Аннотация. Статья посвящена электризации космических аппаратов. В работе рассмотрены два метода расчета картины растекания токов по поверхности космических аппаратов. Первый метод заключается в адаптации наиболее производительного программного обеспечения LTSpice IV для проведения расчетов. Второй метод базируется на разработке редуцированной вычислительной схемы. Приведен сравнительный анализ производительно сти рассматриваемых подходов и программ.
Ключевые слова: программа расчета электрических схем, редукция, система обыкновенных дифференциальных уравнений, сравнительный анализ.
Исследование осуществлено в рамках программы фундаментальных исследований НИУ ВШЭ в 2015 году. На сегодняшний день проблема электризации космических аппаратов (КА) на околоземной орбите стоит достаточно остро. Примерно в 25-30% случаев отказов КА связано именно с этим явлением. Космический корабль в период активного существования пребывает в разреженной плазме. Заряженные частицы, накапливаясь на поверхности КА создают разность потенциалов между различными участками поверхности КА, вследствие чего генерируется электростатический разряд (ЭСР) величиной до 100 А [1]. ЭСР создает наводки во фрагментах бортовой кабельной системы (БКС). Такие электромагнитные помехи способны вывести из строя радиоэлектронную аппаратуру, а значит и сам КА.
Чтобы избежать возможных неполадок в электронике КА необходимо на этапе эскизного проектирования построить картину растекания токов по поверхности КА, рассчитать наводки в БКС, а затем дать рекомендации по оптимальному маршруту прокладки кабелей или их экранированию. С этой целью группой ученых в Московском институте электроники и математики была разработана структурная электрофизическая модель (СЭМ) КА и программное обеспечение (ПО) «Satellite-MIEM» для ее реализации.
СЭМ КА генерируется средствами ПО «Satellite-MIEM» по построенной в 3D Мах полигональной объемной модели. При этом поверхность модели состоит из совокупности прямоугольников, преобразуемых в поверхностную равномерную сетку: совокупность
Программирование методов расчета картины растекания токов
по поверхности космических аппаратов_
связанных узлов. Количество узлов СЭМ крупногабаритного КА равняется (1...2)х105.
Каждую связь (ветвь) можно представить в виде RLC-цепи, образующую, эквивалентную электрическую схему (ЭЭС) поверхности КА. ЭЭС содержит в себе резисторы, катушки индуктивностей, конденсаторы и источник тока. ЭЭС, содержащую вышеперечисленные элементы можно записать в виде системы обыкновенных дифференциальных уравнений. Для расчёта переходных токов в ПО используется программа расчета электрических схем PSpice 9.0 [2].
Построение картины растекания токов является наиболее трудоемким процессом во всей процедуре расчета наводок. Решение данной задачи в ПО PSpice 9.0 для схемы порядка 150000 узлов требует около недели на современных ПК [2], что неприемлемо для предприятий космической отрасли.
Мы предлагаем 2 выхода из сложившейся ситуации. Первый путь заключается во внедрении наиболее производительного ПО. В ходе анализа программ расчета электрических схем [2] было выявлено, что наиболее производительной является ПО LTSpice IV. Поэтому нами было принято решение использовать симулятор LTSpice IV для расчета растекания токов по поверхности КА, внеся определенные изменения в исходный код программы «Satellite-MIEM».
Алгоритм работы программы «Satellite-MIEM» представлен на рисунке 1.
Для подключения нового ядра необходимо было отредактировать исходный код программы так, чтобы формировался новый файл, данные которого будут соответствовать синтаксису новой программы LTSpice IV. Основным отличием в синтаксисе между PSpice 9.0 и LTSpice IV является форма записи резистора, конденсатора и индуктивности, соединений между ними и землей.
Ниже приведена демонстрация двух строк, которые будут записаны в файл для последующего расчета с помощью двух программ на ядре Spice:
RO Gr 0 GR IK R1 (формат PSpice 9.0); R1 N001 0 1000 (формат LTSpice IV)
В обоих примерах описано подключение резистора R1 номиналом 1000 Ом к земле. В соответствии с приведенным примером был модифицирован исходный код программы «Satellite-MIEM» и подключена сторонняя программа LTSpice IV. В результате работы было изменено три основных функции - функция вызова внешнего ядра, функция записи в файл резисторов, соединенных с землей и функция записи в файл резисторов, соединенных между собой; было изменено порядка двухсот строк кода с целью реализовать создание файла с необходимым синтаксисом, понятным LTSpice IV.
Сравнительный анализ производительности программ Р8рюе 9.0 и ЬТ8р1се IV приведен в таблице 1. Эксперименты проводились на ЭВМ с двухъядерным процессором (тактовая частота равняется 2,4 ГГц на каждом ядре) и объемом оперативной памяти, равной 4 ГБ.
Рис.1 Алгоритм работы программы «8а1е11ке-МШМ»
Таблица 1. Сравнительный анализ производительности программ РЗргсе 9.0 и
ЬТЗрюе IV
Количество узлов в ээс Расчетное время Р8рюе 9.0 Расчетное время ЬТБрюе IV
602 3.82 1.56
5652 135.86 14.65
13075 1127.73 68.83
Суть второго пути заключается в разработке нового математического ядра. Идея подхода появилась на основании того, что ток, способный нанести отрицательное воздействие на электронику КА, сосредоточен в локальной области вокруг точки приложения разряда (менее 1% от всей площади поверхности на крупногабаритных КА). Исходя из принципов
Программирование методов расчета картины растекания токов
по поверхности космических аппаратов_
макромоделирования предлагается исключить из математической модели порядка 99% неизвестных, а расчет проводить только в локальной области КА [3,4]. Итак, пусть математическая модель ЭЭС КА сформирована в расширенном однородном координатном базисе в виде системы обыкновенных дифференциальных уравнений
C(Q)^-X{t) + G(Q)X(t) = 7(0, Х(0) = Хо, (1)
dt
где C(Q),G(Q) - (ихп) - матрицы, содержащие в себе номиналы индуктивностей, емкостей и проводимостей эквивалентной электрической
—Т
схемы (ЭЭС), Q ~ <1г) - вектор варьируемых параметров модели,
X(t) _ вект0р ИСкомых фазовых переменных (напряжения во всех узлах
схемы и токи, протекающие через индуктивные элементы), - вектор входных сигналов.
С учетом того, что m<n+nY, исходная задача (1) может быть
записана в блочном виде
Гг г 1 Ml М2 d ~Xi(t) + ~Gn G\2 ~Mt) Ht)
_С21 С22 (Q) dt _X2(t)_ _G2l G22(Q)_ _X2(t)_ J2(t)_
, Хф) = Х0,
(2)
где ^22 (б)'^22 (0 _ (/их/и) _ матрИЦЫ5 Хг(1) _ (/ИХ 1) _ вект0р? содержащий искомые выходные характеристики модели.
Задачу (2) необходимо преобразовать для вычисления только
вектора а так как преобразованная задача будет состоять из ш « п уравнений, трудоемкость вычислений значительно уменьшится.
Для реализации запишем решение системы (2) явным и неявным методом Эйлера:
Явный метод Эйлера задается выражением
Неявный метод Эйлера:
dt
Отсюда получим следующие выражения для исходной системы уравнений (1),
\ сщм) - [С-\ - ошд = Щ), (3)
h
h
[c\ + G]X{tl)-\cX{t1_l) = Y(ti). h h
(4)
В блочном виде система уравнений из двух методов выглядит следующим образом:
"11
"21
'12
"22
Ux 'Сп -hGn г М2 -hGn >1" = h
U2_ г -hG2l с ^22 -hG22_ _V2_ Y20-1)_
(3)
Cn + hGn С2i + hG2l
С12 + hGu
Гс/i" "Qi Г 1 >i" = h Ym
[u2_ r _ 21 с 22 _ /2_ /2(0_
(4)
где векторы и я V - обозначения векторов X и 1,ч соответственно.
В этом случае исходная система дифференциальных уравнений
записана в явной форме и имеет простой вид: •
Х{г) + АЩ) = Ш> Х(0) = Хо. (5)
А явный и неявный методы Эйлера для такой системы в блочном виде записываются следующим образом:
Ux ~Exx~Hx -hA12 >1_ = h Zi"
Ui_ -hA2l E22 ~ _ Z2
Еп + Mi
hA
-21
м2 Е22 ^22.
U1 Vx = h Rx
—
Ui V2 R2
Блоки All, А12,А21 и А22 получаются, исходя из размеров искомой области, как показано на рисунке 2.
All A12
A21 A22
Рис. 2. Преобразование матрицы А
Комбинируя два уравнения из неявного метода и одно из явного, получаем редуцированную формулу для вычисления необходимого подвектора решения. Если возмущающие воздействие в подвекторе У2, то рабочая формула принимает вид:
[В, - кВ2 ]Х2(0 - В, Х2(,-1) = А?2со (6),
А = ^22 +^21(Д1) ~~ Л 2 И —^21(^11) 1 ^11^12 ~Л22
где
Для реализации и проверки значимости редуцированной вычислительной схемы для расчета линейных ЭЭС большой размерности было разработано ПО [5], работу которого можно разбить на следующие этапы:
1. Конвертирование ЭЭС КА из программы ЬТ8рюе в модель, сформированную в РОКБ, необходимую для дальнейших расчетов и преобразований.
Программирование методов расчета картины растекания токов
по поверхности космических аппаратов_
2. Редукция матриц по новой вычислительной схеме для ускорения расчета.
3. Расчет ЭЭС с помощью вычислительной схемы.
Расчет временных затрат. Время анализа ЭЭС, состоящей из 1000 узлов, разными методами при шаге 11=0,001 сек и количестве шагов 1000 приведено в таблице 2, а время построения редуцированной вычислительной схемы приведено в таблице 3. При этом результаты расчета по редуцированной вычислительной схемы отличаются от результатов неявного метода « 1%. Как видно из таблиц 1 и 3, время построения макромодели сравнимо с временем однократного анализа в ПО ЬТ8рюе IV. Таким образом подход, основанный на редукции становится оправданным при многократном использовании макромодели.
Таблица 2. Время расчета ЭЭС явным и неявным методами по сравнению
с редуцированной вычислительной схемой
Численный метод Время расчета (с)
Явный метод Эйлера 11,4
Неявный метод Эйлера 1057,01
Редуцированная вычислительная схема: 3 узла в подвекторе х2 0,063
11 узлов в подвекторе х2 0,069
20 узлов в подвекторе х2 0,073
Таблица 3. Время построения макромодели для редуцированной вычислительной
схемы
Количество узлов в Время построения
подвекторе х2 редуцированной вычислительной схемы (с)
Редуцированная 8,2
вычислительная схема: 3 узла в
подвекторе х2
11 узлов в подвекторе х2 8,4
20 узлов в подвекторе х2 8,9
Практический пример. Так как продемонстрировать работу вычислительной схемы для большого количества узлов не представляется рациональным из-за большого объема данных, для подтверждения результатов работоспособности редуцированной вычислительной схемы приведен пример ЭЭС с небольшим количеством узлов (рисунок 3).
CI CI
I R1 LI ZftS
1» .и"4 !t i#V R2 II
L4 1 N3 u Li M
i < ' ** H
.Vi Mi 1ML12 ' h LI) LI i,.
с с f
j [l
Рис. 3. ЭЭС, содержащая 21 узел
Оценим результаты, полученные после анализа схемы различными вычислительными схемами (рисунок 4). Пусть узел включения к первому узлу подключен источник тока, а искомая локальная область - потенциалы в узлах: 1, 2,4.
X (явный метод Эйлера) (В)
Узел 1:49,9466835702948730002 Узел 2: 0,0016581325264198413 Узел 3: 0,0000000411625455172 Узел 4: 0,0000000000008149999 Узел 5: 0,0000000000000000504 Узел 6: 0,0000000411625454935 Узел 7: 0,0000000683253730828 Узел 8: 0,0000000000022093904 Узел 9: 0,0016581325264198413 Узел 10: 0,0498949354265036122 Узел 11: 0,0000016547350489775 Узел 12: 0,0000016547350489774 Узел 13: 0,0498949354265036122 Узел 14: 0,0000016547621418938 Узел 15: 0,0000000000410362437 Узел 16: 0,0000000000681291605 Узел 17: -0,0000000000000008116 Узел 18: 0,0000000000000022004 Узел 19: 0,0000000000000008113 Узел 20: 0,0000016547621418939 Узел 21: 0,0000000000410376327
X (неявный метод Эйлера) (В)
Узел 1:49,9466735723087110002 Узел 2: 0,0016631308541913806 Узел 3: 0,0000000414949568105 Узел 4: 0,0000000000008274094 Узел 5: 0,0000000000000000513 Узел 6: 0,0000000414949567862 Узел 7: 0,0000000687959991111 Узел 8: 0,0000000000022411892 Узел 9: 0,0016631308541913833 Узел 10: 0,0498950206155334845 Узел 11: 0,0000016597378919406 Узел 12: 0,0000016597378919408 Узел 13: 0,0498950206155334979 Узел 14: 0,0000016597651231802 Узел 15: 0,0000000000413684540 Узел 16: 0,0000000000685996942 Узел 17: -0,0000000000000008244 Узел 18: 0,0000000000000022325 Узел 19: 0,0000000000000008243 Узел 20: 0,0000016597651231807 Узел 21:0,0000000000413698638
Х2(редуцированная вычислительная схема) (В)
Узел 1:49,9466735923236270002 Узел 2: 0,0016631208227460565 Узел 6: 0,0000000414939574496
Рис. 4. Расчет при h=0.001 с и количестве шагов = 1000
При этом, результаты анализа, полученные по редуцированной вычислительной схеме и неявным методом Эйлера отличаются на тысячные доли процента.
Список литературы
[1] Новиков JI.C, Бабкин Г.В., Морозов Е.П., Колосов С.А., Крупников К.К., Милеев
B.Н., Саенко B.C. Комплексная методология определения параметров электростатической зарядки, электрических полей и пробоев на космических аппаратах в условиях их радиационной электризации. Руководство для конструкторов. ЦНИИМАШ, Королев 1995.
[2] Востриков А. В., Абрамешин А. Е. Тестирование коммерческого программного обеспечения для моделирования и анализа эквивалентных электрических схем космических аппаратов// Технологии электромагнитной совместимости. 2012. № 1.
C. 25-28.
Программирование методов расчета картины растекания токов
по поверхности космических аппаратов_
[3] Востриков А. В., Абрамешин А. Е., Борисов Н. И. Расчет наводок в бортовой кабельной сети космических аппаратов с помощью макромоделирования на основе методов Эйлера// Технологии электромагнитной совместимости. 2012. № 1 (40). С. 19-24.
[4] Востриков А. В., Абрамешин А. Е. Вычислительная схема ускоренного метода расчета наводок в бортовой кабельной сети космических аппаратов// Технологии электромагнитной совместимости. 2012. № 3. С. 22-28.
[5] Востриков А. В., Жадов А. Д., Борисов Н. И., Клышинский Э. С. Разработка и анализ редуцированной вычислительной схемы// Системный администратор. 2014. № 12. С. 90-95.
