CC BY
45
Разработка эффективного метода анализа эквивалентных электрических схем космических аппаратов, использующего редукцию и разреженность матриц моделей
Востриков А.В., Борисов Н.И.
МИЭМ, каф. ИТАС. e-mail: sanchs@inbox.ru
О проблеме электризации космических аппаратов (КА) и связанных с ней электростатических разрядах (ЭСР) известно уже около 40 лет [1,2]. ЭСР, которые являются следствием электризации, оказывают негативное влияние на работу бортовой радиоэлектронной аппаратуры (БРЭА) КА. Электромагнитные помехи, создаваемые ЭСР, вызывают сбои в работе различных бортовых систем КА, а интенсивные разрядные токи могут привести к необратимым повреждениям элементов аппаратуры. Основными рецепторами импульсных помех от ЭСР являются фрагменты бортовой кабельной сети (БКС), проложенные по внешней поверхности КА [3].
Для расчета величины наводок в БКС в МИЭМ разработана структурная электрофизическая модель (СЭМ) электризации КА и программное обеспечение (ПО) «Satellite-MIEM» для ее реализации [4,5].
СЭМ КА - эквивалентная электрическая схема (ЭЭС) поверхности КА из сосредоточенных R, L и C элементов. Для построения СЭМ корпус КА с навесными элементами разбивается на такие геометрические тела, как цилиндр, сфера, тор, плоскость, стержень и т.д.
Корпус КА в виде полигональной 3d-модели, состоящей из прямоугольников, можно построить в пакете 3dsMax. Средствами ПО «Satellite-MIEM» полигональная модель преобразуется в поверхностную сетку: совокупность связанных узлов. Каждая связь (ветвь) представляется в виде элементов электрической цепи, в целом образующих ЭЭС поверхности КА.
Зная возможные места возникновения ЭСР с помощью ПО «Satellite-MIEM» можно рассчитать картину растекания токов по поверхности КА. Для расчета величины помехового сигнала во фрагменте БКС средствами ПО прокладывается трасса этого фрагмента по внешней поверхности КА, затем вычисляется величина наводки в кабеле на основе картины растекания токов и коэффициенте трансформации тока протекающего по элементу КА в напряжение наводки. Расчетные данные по уровням помеховых сигналов на входах электронных блоков, полученные на этапе эскизного проектирования КА, включаются в ТЗ на изготовление этих блоков. Таким образом, заказчик получает рекомендации по защите электронных блоков на корпусе КА.
Для расчета картины растекания переходных токов по СЭМ КА ПО «Satellite-MIEM» использует известную программу расчета электрических схем PSpice из пакета OrCAD. Анализ СЭМ, состоящей из 150 000 узлов и более на современных ПК в такой программе невозможен по причине огромной трудоемкости, затрачиваемой ЭВМ на решение системы обыкновенных дифференциальных уравнений (ОДУ) методом Рунге-Кутта. Исследования показали, что время, затраченное на решение системы 2000 и более ОДУ в ПО PSpice, неприемлемо для расчетов [6]. Возникла необходимость проводить анализ СЭМ КА с адекватными временными затратами.
Авторами предлагаются 2 выхода из создавшейся ситуации.
Особенностью методов является проведение расчетов переходных токов лишь в локальной зоне схемы КА. Возможность использование таких подходов оправдано. Проведенные эксперименты в программе «8а1еШ1е-М1ЕМ» показывают, что область, в которой значения переходных токов в ветвях не больше величины в 1-2 % от величины ЭСР, представляет собой 400 узлов вокруг места возникновения ЭСР. А значения переходных токов в дальних ветвях, не превышающие величины в 1-2 % от амплитуды тока в точке приложения ЭСР, считаются незначительными в связи с тем, что электромагнитные помехи от таких переходных токов оказывают слишком слабое влияние на БРЭА КА.
В статье [6] был предложен метод приближенного ускоренного расчета растекания токов от электростатического разряда по поверхности космического аппарата с помощью программы Р8р1ее. Суть метода состоит в построении ограниченной области расчета растекания токов. Расчет переходных токов будет происходить только в этой, обозначенной пользователем зоне СЭМ, не затрагивая ее остальную часть.
В работе [7] описан другой подход, исключающий недостаток приближенного метода ускоренного расчета, заключающийся в использовании стороннего ПО. Идея подхода заключается в построении макромодели на основе явного и неявного методов Эйлера. По макромодели будут найдены необходимые неизвестные (потенциалы в узлах или токи в ветвях). Основные шаги вывода новой вычислительной схемы следующие:
1) Формулировка модели схемы в РОКБ особым образом:
(С„ 0 0 1 Ь, 0 0
"11 0
0
'22 0
) ( О О12 О,3 ' ~ X!«) " )
X2 ( ) + О21 0 О23 X2 ( ) = 0 , (1)
) V О31 О32 О33 у _ Xз(t) _ 0
где все матрицы симметричные, подматрицы С11, Ь22, О33 - диагональные, при этом:
- X1 (^) - подвектор напряжений в узлах, к которым подсоединены конденсаторы;
- Х2 (^) - подвектор токов, проходящих через индуктивные элементы;
- X3 (^) - подвектор напряжений в узлах, к которым не подсоединены конденсаторы.
2) Исключение из модели подвектора ). Поскольку матрица О33 диагональная, а остальные подматрицы сильно разрежены, то получим следующую модель с разреженными матрицами (подробный вывод приведен в [7]):
(С 0 >1 0 4
2 У
Х() 4(0
+
V 21
аl3гзза31 О12 -О г О 1 13 33 32 ~Xl(t)'
а2зг3за3l О23 - а2згззаз2 у _X2(í)_ 0
(2)
где т33 - сопротивление резисторов в схеме.
С учетом того, что С11 и Ь22 - диагональные матрицы, модель может быть записана в виде системы дифференциальных уравнений, заданной в явной форме:
X (0) = X 0, (3)
^X (?) + АХ V) = У (0,
где А =
(о - О г О О - О г О 1( С 0 1
^11 1-г13/Э3^731 12 13 33 32 и
V О21
"О23Г33°31
а
23
"а23Г33а32
0 ь
'22 У
Полученная матрица А может получиться как плотной в случае небольших СЭМ КА или схем другой структуры, так и разреженной в случае большой эквивалентной электрической схемы.
3) Запись явного и неявного вычислительных процессов в блочной форме:
(4)
(5)
Ц1" _ 'V1" " А11 А12 'V1" Ь +
Ц 2 _ V 2 _ _ А21 А22 _ V 2 _
Еп + Ып
ЬА
ЬА
12
E22 + ЬА22,
U1 Ц 2
V1
V 2
У1 Ь,
У 2 _
У1
+
У 2
h,
где [ц 1 Ц2 ] = XTI+l, [V1 V2 ] = X] .
Будем считать, что подвекторы Ц 2 и V2 содержат по т << п искомых коэффициентов решения. Необходимо на базе выражений (4) и (5) сформулировать новую вычислительную схему, из которой будут исключены подвекторы Ц1 и V1. Для этого запишем (4) и (5) в виде следующих подсистем уравнений
Ц1 = (Е11 - hA11)V 1 - ЬА12 V 2 + hY 1, Ц2 = -ЬА21У 1 + (Е22 - ЬА22)У2 + ЬУ2, (Е11 + ЬА11)Ц 1 + ЬА12Ц 2 = V1 + ЬУ1,
(4.1)
(4.2)
(5.1)
(5.2)
ЬА21Ц1 + (Е22 + ЬА22)Ц 2 = V 2 + ЬУ 2.
4) Далее исключим из модели подвекторы Ц1 вычислительную схему (6), построенную с использованием макромоделирования. Подробный вывод формулы рассмотрен в [7].
и V!,
получим новую идеи
[Е22 + А21(А^Г142]Ц2 = [А^ГЕ + ЪАп)А12 + (Е22 -ЪА22)]У2 + Ь[-А21(А„)-1 АиУ 1 + У2]
(6)
После ряда преобразований получим вычислительную схему:
(X 2)г+1 = [СЬ + С2](Х 2)г + ЬСз(У 2), ] (7)
В итоге была построена макромодель и получена новая вычислительная схема, в которой вычисляются только т << п коэффициентов вектора X(^).
В случае больших схем специфика задачи состоит в том, что необходимо решать систему линейных алгебраических уравнений (СЛАУ) с разреженной блочной матрицей А высокого порядка. Наибольшей проблемой в (6) является обращение числовой подматрицы А11 высокого порядка. Трудоемкость нахождения обратной
матрицы стандартным QR-разложением равна п3 ВМО. Таким образом, если, например, А11 имеет размерность (150000*150000), то для ее обращения потребуется
3375 *1012 вещественных мультипликативных операций (ВМО). Необходимо принять меры по сокращению временных затрат на обращение матрицы. Поскольку матрица А разрежена, предлагается преобразовать к треугольному виду с окаймлением с помощью метода определяющих величин [8] либо подматрицу А11 либо матрицу А.
Требования к формированию блочной матрицы треугольного вида с окаймлением:
1) Размер окаймления должен быть минимальным. Чем меньше размер окаймления, тем быстрее будет решаться СЛАУ.
2) Элементы, стоящие на диагонали должны по модулю быть как можно больше. Дело в том, что в процессе решения СЛАУ приходится обращать треугольную матрицу. Чем больше значения по модулю на диагонали, тем больше определитель матрицы, что влечёт за собой высокую точность обращения матрицы.
Алгоритм преобразования разреженной блочной матрицей А высокого порядка к треугольному виду с окаймлением с помощью метода определяющих величин сводится к перестановке строк и столбцов преобразуемой матрицы и приведен в [8]:
Шаг 1. Присвоить индексу перебора строк г значение единицы (г = 1).
Шаг 2. Определить строки матрицы А, которые связаны с минимальным числом переменных. Пусть минимум этого числа равен к.
Шаг 3. Если к = 1, то перейти к шагу 5.
Шаг 4. Выбрать строку, которая связана с к переменными и сумма ненулевых элементов в столбцах которой максимальна.
Шаг 5. Пусть выбрана ]-я строка. Переставить 1-ю и ]-ю строки.
Шаг 6. Все выбранные к-е переменные исключить из дальнейшего рассмотрения (так как они уже определены). Проверить, появились ли строки, которые не связаны ни с одной из оставшихся переменных. Если такие строки есть, то исключить их из дальнейшего рассмотрения, при этом г = г +1.
Шаг 7. Если перебраны не все строки, тогда изменить индекс г на единицу (г = г +1) и возвратиться к шагу 2.
Замечание. Применение метода рассматривалось ранее в [9] для вычисления собственных значений полиномиальной матрицы. В данной работе метод определяющих величин используется для решения числовой системы ОДУ большой размерности.
После работы алгоритма матрица А разбивается на блоки. Трудоемкость Т процесса преобразования квадратной разреженной (п*п) - матрицы А к треугольной матрице с окаймлением в худшем случае, когда нужно переставить все строки и все столбцы:
Т »(п -1)2 ВМО.
Трудоемкость обращения треугольной матрицы равна 1 п3 ВМО.
Таким образом, процесс нахождения обратной матрицы сводится к (п -1)2 +1 п3 ВМО.
Но в вычислительной схеме (6) трижды присутствует произведение А21( А112)-1.
Общая трудоемкость вычисления произведения без преобразования
Т = 2п3 + п2т ВМО.
Для снижения трудоемкости расчета по вычислительной схеме (6) предлагается не вычислять обратную матрицу (А112)-1, а найти произведение А21(А112)-1, решив СЛАУ:
(А21(АП2)-1)* Ап2 = А21 (8)
Если подматрица А11 имеет размерность (п*п), а у подматрицы А21 размерность (т*п), то трудоемкость нахождения неизвестной величины:
Т = (п -1)2 + дт + д + notdiag ВМО,
где д - количество ненулевых элементов подматрицы А11, по1&а§ - количество недиагональных ненулевых элементов подматрицы А11. При п = 150000, т = 400, д = 1000000 расчет от нахождения обратной матрицы будет отличаться в 295129 раз. Таким образом, речь идет об экономии машинного времени на 5 порядков.
Ранее указывалось о двух подходах к решению задачи:
- В первом случае предлагалось преобразовать подматрицу Л11 к треугольному виду с окаймлением;
- Во втором предлагалось провести преобразования с матрицей А.
Характеристика первого подхода.
После преобразования матрицы А необходимо получить блоки: Л11 - не дефектная нижне-треугольная матрица с окаймлением порядка (п*п), Л12 матрица (п*т), Л21 матрица (т*п), Л22 матрица (т*т) (рис. 1). При этом необходимо, чтобы подматрица Л22 содержала в себе такие номиналы элементов, узлы и ветви которых будут присутствовать в макромодели.
* * * * * * * * * * *
* * * * * * * *
* * * * * * * * * * 0 * * *
* * * * * * * * к
* * * А * * * * \ А А11 * А12
* * * * * * * * V * * * * * *
* * * * * * * * * * * * * * *
* * * * * * * * * * * * * * * *
* * * * * * * * * * * * * * * *
* * А21 * * * А22
Рис. 1. Преобразование подматрицы Л11 к треугольному виду с окаймлением
Трудоемкость нахождения Л21( Л112)-1 Т сводится примерно к (п -1)2 + цш + ш2п + ц + notdiag ВМО.
Преимущество метода:
- нет необходимости в действиях после работы вычислительной схемы
Недостаток:
- представляется трудоемкой задача нахождения произведения Л21( Л112)-1.
Характеристика второго подхода.
После преобразования матрицы А получим блоки: Л11 - не дефектная нижнетреугольная матрица порядка (п*п), Л12 матрица (п*т), Л21 матрица (т*п), Л22 матрица (т*т) (рис. 2).
* * * * * * * * * *
* * * * * * * *
* * * * * * * * * * 0 * *
* * * * * * * * к
* * * А * * * * \ У А11 А12
* * * * * * * * V * * * * * *
* * * * * * * * * * * * *ч * *
* * * * * * * * * * * * * * * *
* * * * * * * * * * * * * * * *
* * А21 * * * А22
Рис. 2. Преобразование матрицы А к треугольному виду с окаймлением Покажем пример, реализованный на ЭВМ. Рассмотрим плоскость, состоящую из двух ячеек. Номиналы всех элементов в ЭЭС равны 1 для простоты. Матрица А до и после преобразования показана на рис. 3.
-1 0 1 0 0 0 0 0 0 1 0 0 0 -1 0 0 0 0 0 0 0 0 0 0 0 -1
1 2 -1 0 1 0 -1 1 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 -1 0
-1 0 0 0 0 0 0 0 0 0 0 -1 0 -1 1 1 0 0 0 0 0 0 0 0 0 0
0 0 0 1 1 0 0 0 1 0 0 1 0 0 -1 0 1 0 0 0 0 0 0 -1 0 0
0 0 0 0 1 0 0 -1 -1 0 -1 0 0 -1 0 0 0 -1 0 0 0 0 1 0 0 0
0 0 0 1 0 2 1 0 0 -1 1 0 0 к 0 -1 0 0 -1 1 0 0 0 0 0 0 0
\
0 0 1 0 0 -1 0 0 0 0 0 0 0 \ 0 0 1 0 0 1 -1 0 0 0 0 0 1
0 0 0 0 1 -1 0 0 1 1 0 0 0 V 0 0 0 0 0 -1 -1 -1 0 1 0 0 0
0 0 0 -1 1 0 0 0 0 1 0 1 0 0 0 0 -1 -1 1 0 0 0 -1 0 0 0
0 0 0 -1 -1 0 0 0 0 0 1 0 -1 0 0 0 0 1 1 1 0 0 0 0 0 1
0 -1 -1 0 0 0 0 1 0 0 0 0 0 0 0 1 0 1 1 0 0 0 0 0 -1 0
-1 0 0 0 0 0 0 0 -1 0 1 0 0 0 0 -1 0 0 0 1 0 1 1 0 2 0
0 0 -1 0 -1 0 0 0 -1 0 0 0 0 1 -1 0 1 0 1 0 0 -1 0 2 0 0
Рис. 3. Преобразование матрицы модели плоскости ЭЭС, состоящей из двух ячеек к треугольному виду с окаймлением Преимущество метода:
- на вычисление А21(А112)-1 потребуется меньшее количество машинного времени.
Недостаток:
- дальнейшие вычисления после работы вычислительной схемы.
Пояснение недостатка. После преобразования к треугольному виду с
и А
окаймлением с помощью метода определяющих величин матрицы А в окаймление могут попасть как номиналы элементов (Х2)г.+1, узлы и ветви которых присутствуют в макромодели так и элементы не присутствующие в макромодели.
Рис. 4. Анализ модели после преобразования матрицы А
Проделав по вычислительной схеме необходимое число шагов, будет вычислен подвектор (Х2);+1. Потенциалы в узлах необходимой локальной области вычисляются исходя из следующих соображений (рис. 4):
АХ = У ,
или
(А1 + А21) X + (А12 + А22) Х2 = У2,
или
Ап х+А12 х2 = у
или
X! = - ^12 X2) В силу того, что на практике подвектор Y1 нулевой, получим:
Xi = (-V A12) X2
трудоемкость
Произведение Л11 Л12 вычисляется аналогично (8), тогда вычисления неизвестного подвектора:
Т » цш + шп ВМО.
Общая трудоемкость вычисления Л21( Л112)-1 без преобразования
Т = 2п3 + п2ш ВМО. Общая трудоемкость после преобразования
Т = (п -1)2 + 2цш + шп + ц + notdiag ВМО,
или
Т » п2 + 2цш + шп + ц + notdiag ВМО.
При п = 150000, п = 400, q = 1000000 на вычисление
преобразований уйдет 6759*1012 ВМО, а после преобразования 2336185*104 ВМО. Экономия памяти составит примерно 5 порядков.
Выводы
1. Было показано, что в случае больших схем необходимо решать СЛАУ с разреженной блочной матрицей А высокого порядка. Учитывая данную специфику, было предложено преобразовать эту матрицу к треугольному виду с окаймлением с помощью известного метода определяющих величин.
2. Рассмотрены случаи преобразования к специальному виду матрицы А. Выделены преимущества и недостатки методов.
A21( A11 )
без
Список литературы
1. Garrett H.B. The Charging of Spacecraft Surfaces / Review of Geophysics and Space Physics. 1981. V. 19. № 4. P. 577-616.
2. Тютнев А.П., Саенко В.С., Пожидаев Е.Д., Сасов А.М. Моделирование электрических явлений, сопровождающих воздействие низкоэнергетической компоненты космических ионизирующих излучений // Конструирование научной космической аппаратуры. М.: Наука, 1982. С.66-80.
3. Новиков Л.С., Бабкин Г.В., Морозов Е.П., Колосов С. А., Крупников К.К., Милеев В.Н., Саенко В.С. Комплексная методология определения параметров электростатической зарядки, электрических полей и пробоев на космических аппаратах в условиях радиационной электризации. Руководство для конструкторов. -ЦНИИмаш, 1995. - 160 с.
4. Марченков К.В., Соколов А.Б., Саенко В.С., Пожидаев Е.Д. Новое поколение программного обеспечения «Satellite-MIEM» для расчета наводок во фрагментах бортовой кабельной сети, проложенных по внешней поверхности космических аппаратов // Технологии электромагнитной совместимости. - Москва, изд-во ООО «Издательский Дом «Технологии». - 2008. - № 1(24), С. 39-44.
5. Марченков К.В., Дорофеев А.Н., Востриков А.В., Саенко В.С. Новое поколение программного обеспечения «Satellite-MIEM» для расчета наводок во фрагментах бортовой кабельной сети, проложенных по внешней поверхности
космических аппаратов. Труды XVII Международного совещания «Радиационная физика твердого тела», 30 июня - 5 июля 2007 г. Севастополь. - С. 421-425.
6. Востриков А.В. Приближенный метод расчета растекания токов по элементам конструкции космического аппарата при электростатических разрядах. Технологии ЭМС №2(33). Москва 2010, с. 75-79.
7. Востриков А.В., Борисов Н.И. Новый алгоритм построения макромоделей на основе методов Эйлера. Труды XXI Международного совещания «Радиационная физика твердого тела», 22 августа - 27 августа 2011 г. Севастополь, с. 283-291.
8. Петренко А.И., Власов А.И., Тимченко А.П. Табличные методы моделирования электронных схем на ЭЦВМ. - Киев: Вища школа, 1977. - 192 с.
9. Баскаков А.Е. Алгоритм приведения разреженной матрицы к треугольной форме с окаймлением в рамках макромоделирования - Сб. трудов научно-практического семинара «Новые информационные технологии в автоматизированных системах», М., 2009.
