Программирование методов расчета картины растекания токов по поверхности космических аппаратов

Программирование методов расчета картины растекания токов по поверхности космических аппаратов

Борисов Н.И., Востриков А.В., Жадов А.Д., Колтунова В.Э.

Национальный исследовательский университет Высшая школа экономики 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

В обоих примерах описано подключение резистора 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

dt

где C(Q),G(Q) - (ихп) - матрицы, содержащие в себе номиналы индуктивностей, емкостей и проводимостей эквивалентной электрической

—Т

схемы (ЭЭС), Q ~ <1г) - вектор варьируемых параметров модели,

X(t) _ вект0р ИСкомых фазовых переменных (напряжения во всех узлах

схемы и токи, протекающие через индуктивные элементы), - вектор входных сигналов.

С учетом того, что m<n+nY, исходная задача (1) может быть

записана в блочном виде

Гг г 1 Ml М2 d ~Xi

_С21 С22 (Q) dt _X2(t)_ _G2l G22(Q)_ _X2(t)_ J2(t)_

, Хф) = Х0,

где ^22 (б)'^22 (0 _ (/их/и) _ матрИЦЫ5 Хг(1) _ (/ИХ 1) _ вект0р? содержащий искомые выходные характеристики модели.

Задачу (2) необходимо преобразовать для вычисления только

вектора а так как преобразованная задача будет состоять из ш « п уравнений, трудоемкость вычислений значительно уменьшится.

Для реализации запишем решение системы (2) явным и неявным методом Эйлера:

Явный метод Эйлера задается выражением

Неявный метод Эйлера:

dt

Отсюда получим следующие выражения для исходной системы уравнений (1),

\ сщм) - [С-\ - ошд = Щ),

h

h

[c\ + G]X{tl)-\cX{t1_l) = Y(ti). h h

В блочном виде система уравнений из двух методов выглядит следующим образом:

"11

"21

'12

"22

Ux 'Сп -hGn г М2 -hGn >1" = h

U2_ г -hG2l с ^22 -hG22_ _V2_ Y20-1)_

Cn + hGn С2i + hG2l

С12 + hGu

Гс/i" "Qi Г 1 >i" = h Ym

[u2_ r _ 21 с 22 _ /2_ /2(0_

где векторы и я V - обозначения векторов X и 1,ч соответственно.

В этом случае исходная система дифференциальных уравнений

записана в явной форме и имеет простой вид: •

Х{г) + АЩ) = Ш> Х

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

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.