CC BY
103
Разработка алгоритмов, основанных на редукции математических моделей электрических схем для расчета растекания токов по элементам конструкции космических аппаратов при электростатических разрядах
Борисов Н.И., Востриков А.В.
МИЭМ, каф. ИТАС Контактная информация: тел.: 8-926-566-35-50, e-mail: [email protected]
1. Введение
Одним из факторов, ограничивающих надежную и длительную эксплуатацию космических летательных аппаратов (КЛА) на высокоэллиптических и геостационарной орбитах, является электризация и связанные с ней электростатические разряды (ЭСР) на борту космического аппарата. ЭСР вызывают сбои в работе бортовой радиоэлектронной аппаратуры [1].
Попытки полного исключения возможности возникновения ЭСР путем подбора материалов внешней поверхности космического аппарата (КА) или активной защиты КА до настоящего времени успехом не увенчались. Удается лишь снизить частоту и мощность ЭСР, но не исключить их полностью [2]. Поэтому необходимо принимать дополнительные меры для безотказной работы электроники КА, при воздействии на нее ЭСР.
Одним из важнейших способов исключения возможных негативных последствий для КА, являющихся результатом воздействия электростатического разряда, представляется моделирование картины растекания токов по поверхности КА и расчет величин возникающих помеховых сигналов во фрагментах БКС, проложенных по внешней поверхности КА, еще на этапе проектирования этого аппарата [1]. С этой целью в МИЭМ была разработана структурная электрофизическая модель (СЭМ) электризации КА и программное обеспечение (ПО) «Satellite-MIEM» для ее реализации.
Структурная электрофизическая модель (СЭМ) электризации космических аппаратов используется для получения картины растекания токов по конструкции КА от электростатических разрядов на поверхности этих КА. Эта картина
растекания токов служит в дальнейшем исходной информацией для расчета электромагнитных наводок во фрагментах бортовой кабельной сети, проложенных по внешней поверхности КА.
Принципы построения СЭМ КА базируются на представлении космического аппарата в виде эквивалентной электрической схемы из сосредоточенных Я, Ь и С элементов. Такой подход, в частности, широко и плодотворно применяется при исследовании различных тепловых и электрических процессов. Анализ СЭМ позволяет произвести расчет переходных токов в корпусе КА и указать участки поверхности, вблизи которых нежелательно размещение фрагментов бортовой кабельной сети (БКС). Кроме того, на этапе эскизного проектирования КА разработчики электронных блоков получают параметры импульсных помех в БКС на входах этих блоков.
В качестве входных данных ПО «8а1еШ1е-М1ЕМ» принимает трехмерную (3Б) модель КА, состоящую из совокупности элементарных фигур - треугольников (Такую модель можно подготовить, например, в пакете 3д8Мах). Средствами программы набор треугольников преобразуется в поверхностную сетку (совокупность связанных узлов). Каждая связь представляется в виде элементов электрической цепи, в целом образующих эквивалентную электрическую схему поверхности КА (рис.1)
[3].
Рис. 1. Преобразование совокупности треугольников в СЭМ Зная возможные места возникновения ЭСР, можно рассчитать картину растекания переходных токов по корпусу и навесным элементам КА. При этом ЭСР представляются в виде импульсного источника тока с заранее заданными характеристиками, соответствующими параметрами разряда. Анализ СЭМ позволяет произвести расчет переходных токов в корпусе КА и указать участки поверхности,
вблизи которых нежелательно размещение фрагментов БКС. Кроме того, на этапе эскизного проектирования КА разработчики электронных блоков получают параметры импульсных помех в БКС на входах этих блоков.
Для получения картины растекания токов по поверхности КА с высокой степенью точности необходимо решать систему из 15GGGG обыкновенных дифференциальных уравнений. Для расчета переходных токов синтезированной СЭМ ПО «Satellite-MIEM» использует программу расчета электрических схем SPICE. SPICE-подобные программы (PSpice, HSpice, LTspice, MicroCap, Circuit Maker, Dr. Spice, View Spice и др.) основаны на машинном составлении системы обыкновенных дифференциальных уравнений электрической цепи и их решении без применения упрощающих процедур. В программах используются численные методы Рунге-Кутта или метод Гира для интегрирования системы дифференциальных уравнений.
Результат растекания токов на больших эквивалентных электрических схемах (ЭЭС) (количество узлов схемы равно 15G GGG) на современных персональных ЭВМ либо не может быть получен (случай RLC-схемы), либо расчет занимает очень много машинного времени. Исследования показали, что время, затраченное на решение системы 2GGG и более обыкновенных дифференциальных уравнений в ПО PSpice, неприемлемо для расчетов [4]. Таким образом, применение программного обеспечения PSpice не отвечает требованиям решения больших задач на персональном компьютере. Поэтому возникает необходимость сокращения времени расчета растекания токов по поверхности КА.
В данной работе для решения типовых задач предлагается ввести новую вычислительную схему, основанную на явном и на неявном методе Эйлера.
2. Постановка задачи
2.1. Модель схемы, сформированная в расширенном однородном координатном базисе (РОКБ), может быть записана в виде системы линейных обыкновенных дифференциальных уравнений
C—X (t) + gX (t) = Y (t), X (G) = X g, (2.1)
dt
где С, G - числовые (п х п) - матрицы порядка, X(/) - вектор искомых фазовых переменных (напряжений во всех узлах схемы и токов, протекающих через индуктивные элементы), У (/) - вектор входных сигналов.
Отметим специфику модели схемы.
1). Матрица С является диагональной и вырожденной. Диагональ состоит из 2-х групп коэффициентов и коэффициенты внутри каждой группы одинаковы.
2). Матрица G - невырожденная, симметричная и разреженная матрица.
3). Вектор У(/) содержит только один ненулевой коэффициент вида у/, где /
может быть любым числом от 1 до п .
Сформируем модель схемы в РОКБ следующим образом:
( С М1 0 0 ^ Д(/) (С11 С12 ^13 ^ " Х&)" "Ц/)"
0 ^22 0 Х2(/) + С21 0 С23 Х2(/ ) = 0
V 0 0 0, ^^3(/) V С31 С32 С33 _ Х3(/) _ 0
Отметим, что матрицы симметричные, подматрицы Сп , ^2, - диагональные, при
этом:
- Хх(/) - подвектор напряжений в узлах, к которым подсоединены конденсаторы;
- Х2(/) - подвектор токов, проходящих через индуктивные элементы;
- X 3(/) - подвектор напряжений в узлах, к которым не подсоединены конденсаторы. Модель (2.2) можно записать в виде:
С11X1 (/)+оп X (/)+о12 Х2 (/)+о1Ъ Х3 (/) = У1(/), (2.3)
^22 Х2 (^) + С21Х1 (^ ) + С23 Х3 (^) = 0 , (2 4)
С31Х1(/) + С32Х2 (/) + С33Х3(/) = 0 , (2.5)
Исключим из модели подвектор Х3(/). Из (2.5) получим:
Хъ(/) = -с^Ч^ Х,(/) + С32 Х2(/)) , (2.6)
Подстановка Х3(/) из (2.6) в (2.3) и (2.4) приводит к следующим выражениям:
Г C о ^ 1 •I&T 1 Г + V
10 ^22 у _Х2«
G13G31r G2 —1ВД2 r ' " Х^) w
G23G31r G23 —G23G32ry _ Xl(‘)_ 0
C11 X1 (t) + (G11 — G13G33 G31)X1 (t) + (G12 - G13G—3 G32 )X2 (t) = Y1(t), (2.7)
L22X2 (t) + (G21 - G23G33G31 )X1 (t) + (G23 - G23G33G32 )X2(t) = 0, (2.8)
Поскольку матрица G33 диагональная, а остальные подматрицы сильно разрежены, то получим следующую модель с разреженными матрицами:
(2-9),
3^32' y|_X2()J L0 J
где г - сопротивление резисторов в схеме.
С учетом того, что C11 и L22 - диагональные матрицы, модель может быть записана в виде системы дифференциальных уравнений, заданной в явной форме:
X (t) + AX (t) = Y (t), (2.10)
Необходимо вычислить с минимальными временными затратами решение системы уравнений в момент времени t , т. е. вычислить числовой вектор X(t ).
3. Обсуждение способов решения задачи
*
3.1. Величина t больше длительности переходного процесса системы уравнений.
*
Если заранее известно, что величина t превышает длительность переходного процесса, то решение системы (2.1) совпадает с ее частным решением вида
R(t) = R1t + Ro. Это решение может быть вычислено тривиальным образом путем решения двух систем линейных алгебраических уравнений (СЛАУ) с одинаковыми матрицами
GR1 = Y, (3.1.1)
GR о = —CR1. (3.1.2)
Таким образом, искомое решение будет иметь вид X(t ) = R1t + Rо. В данном случае можно обсуждать лишь выбор метода решения систем (3.1.1) и (3.1.2), использующего симметрию, разреженность и специальную структуру (скорее всего, ленточную) матрицы G .
Замечание. Желательно знать точную структуру и норму матрицы G для оценки сходимости и скорости сходимости итерационных методов решения СЛАУ.
*
3.2. Величина г меньше длительности переходного процесса системы уравнений.
В этом случае необходимо решать систему (2.1). Возможно использование различных методов решения задачи, численных, численно-аналитических, основанных на макромоделировании и т. д.
4. Оценка эффективности использования численных методов
Рассмотрим, каким образом может вычисляться решение системы (2.10), состоящей только из дифференциальных уравнений.
4.1. Оценка снижения трудоемкости процесса решения за счет редукции модели (исключения некоторого количества уравнений).
Будем считать, что исходная задача имеет следующий вид:
—X (г) + АХ (г) = у (г), X (0) = X о. (4.1.1)
При использовании численных методов процесс решения заключается в
*
вычислении вектора решения в моменты времени г0,г1,г2,...,г , т. е.
X(го),Х(*1),...,X(f). При этом решение в точке вычисляется по решению в точке -1 (одношаговый метод) или по решениям в нескольких предыдущих точках (многошаговый метод). Решения во всех точках (кроме первой го = 0, для которой
заданы начальные условия (4.1.1) X(0) = X0) вычисляются с погрешностью [5].
Если погрешность решения от шага к шагу увеличивается (поскольку текущее решение определяется по предыдущим решениям, вычисленным с погрешностью) и начинает неограниченно возрастать, то метод будет численно неустойчивым.
Рассмотрим простейший одношаговый метод Эйлера отыскания решения задачи (4.1.1).
Явный метод Эйлера решения системы уравнений
-&Х(г) = )], Х(0) = X 0, (4.1.2)
&
может быть получен на основе разложения решения в усеченный ряд Тейлора в окрестности к какой-либо точки
— — & — о
X (г. + к) = X (г.) + - X (г. Щ + о(%). (4.1.3)
аг
Подставив в (4.1.3) вместо производной вектора решения его выражение из (4.1.2), получим формулу явного метода Эйлера
X(tм) — X(tl) + Щ, (4.1.4)
где г.+1 = + к, к - величина текущего шага по оси времени, локальная
погрешность интегрирования (погрешность на шаге без учета погрешности решения на предыдущих шагах) равна остаточному члену ряда Тейлора и примерно равна
второй производной решения, умноженной на кI. Ясно, что чем меньше величина Ц, тем меньше локальная погрешность интегрирования (и тем больше шагов по оси
*
времени необходимо сделать для достижения точки г ). Более простая форма записи явного метода Эйлера имеет вид
XI+1 - XI + Е(XI )к. (4.1.5)
С учетом конкретного вида системы уравнений (4.1.1) получим следующую формулу метода Эйлера
XI+1 — XI + (—AXI + УI )к. (4.1.6)
Из [5] известно, что для обеспечения численной устойчивости явного метода Эйлера величина шага не должна превышать абсолютной величины минимального собственного значения матрицы А. В практических задачах именно это условие
*
может приводить к тому, что для достижения точки г из точки г0 потребуется
2 4
сделать 102 —104 шагов по формуле (4.1.6).
Трудоемкость каждого шага примерно равна количеству вещественных мультипликативных операций (операций умножения и деления) (ВМО), необходимых
для умножения матрицы А на текущий вектор решения X. в (4.1.6). Если (п х п) -матрица А плотная, то трудоемкость шага примерно равна п (ВМО), если
разреженная, то трудоемкость равна количеству ненулевых коэффициентов матрицы
А.
Рассмотрим неявный метод Эйлера, сохраняющий численную устойчивость при любой величине шага. Он выводится также из разложения решения в усеченный ряд Тейлора, но отличающегося от (4.1.3)
— — а — о
^ — к ) = ^ ) — -X (г.) к + о(к). (4.1.7)
аг
Аналогично (3.1.4) получим выражение
X (г.—о — X (г.) — Е [ X (г. )]к,
или
X (г.) — X (г.—,) — Е [ X (г. )]к,
или
X.+1 — X. + F(Xi+\)hi. (4.1.8)
Отличие неявного метода Эйлера (4.1.8) от (4.1.5) заключается в том, что
искомое решениеXi+l находится в правой и левой частях выражения (4.1.8). Его не
удается явно выразить через X.. Для вычисления Xi+l (в отличие от (4.1.5)) необходимо решить систему алгебраических уравнений (4.1.8).
С учетом конкретного вида системы уравнений (4.1.1) получим следующую формулу метода Эйлера
X.+1 — X. + (—AX.+1 + У. )к. (4.1.9)
или
(Е + к. А) X.+1 — X. + У .к.. (4.1.10)
Как следует из (4.1.10), на каждом шаге необходимо решать СЛАУ. При этом
величина шага ограничивается только допустимой локальной погрешностью интегрирования. Шаг может изменяться (обычно, увеличиваться) в зависимости от поведения решения [5].
Если (п х п) - матрица А плотная, то трудоемкость шага примерно равна п3
(ВМО) при каждом изменении величины шага к. Если величина шага не изменялась,
то трудоемкость шага примерно равна п2 (ВМО). Если матрица А разреженная, то
2
трудоемкость шага примерно в п раз меньше чем в случае плотной матрицы.
Рассмотрим возможность снижения трудоемкости вычислений путем предварительного преобразования задач (4.1.6) и (4.1.10) в некоторую новую задачу, количество уравнений в которой будет равно т << п . Особенностью метода является проведение расчетов переходных токов лишь в локальной зоне схемы КА. Проведенные эксперименты в программе «8а1еШ1е-М1ЕМ» показывают, что область, в которой значения переходных токов в ветвях не больше величины в 1-2 % от величины ЭСР, представляет собой 400 узлов вокруг места возникновения ЭСР [4]. А значения переходных токов в ветвях, не превышающие величины в 1-2 % от величины ЭСР, считаются незначительными в связи с тем, что электромагнитные помехи от таких переходных токов оказывают слишком слабое влияние на бортовую радиоэлектронную аппаратуру КА. Поэтому для ускоренного расчета растекания токов от электростатического разряда по поверхности космического аппарата рационально использовать предложенный в настоящей работе метод. Для этого запишем оба вида вычислительных процессов в блочной форме:
У г
У 2
и 1" "V1" " А11 А12 "V1" к +
= —
и 2 _ V 2 _ _ А21 А22 _ V 2 _
к.
(4.1.11)
Е11 + кА11
кА
21
кА12 Е22 + кА22.
и 1
и 2
V1 V 2
Т
+
У1 У 2
к.
(4.1.12)
где [и 1 и 2 ] = X. +1, [V1 V 2 ] = X. . Специфика состоит в том, что матрица А
разрежена и имеет высокий порядок. В связи с этим предлагается преобразовать матрицу А к треугольному виду с окаймлением с помощью метода определяющих
величин. Здесь А11 является не дефектной нижне-треугольной матрицей с не равными нулю диагональными элементами.
Будем считать, что подвекторы и2 и V 2 содержат по т << п искомых коэффициентов решения. Необходимо на базе выражений (4.1.11) и (4.1.12) сформулировать новую вычислительную схему, из которой будут исключены
подвекторы и 1 и V1.
Для этого запишем (4.1.11) и (4.1.12) в виде следующих подсистем уравнений
и 1 = (Ец — кАц) V1 — кАи V 2 + кУ 1, (4.1.11.1)
и 2 = — кА21 V1 + (Е22 — кА22) V 2 + кУ 2, (4.1.11.2)
(Ец + кАп)и 1 + кАпи 2 = V1 + кУ 1, (4.1.12.1)
кА21 и 1 + (Е22 + кА22 )и 2 = V 2 + кУ 2. (4.1.12.2)
Выражения (4.1.11.1), (4.1.11.2), (4.1.12.1), (4.1.12.2) будем рассматривать как
аналог системы из 4-х уравнений с четырьмя неизвестными и 1,и2, V1 ,V2, опирающейся, вообще говоря, на вычислительные процессы разных типов (явного и неявного).
Подставим (4.1.11.1) в (4.1.12.1):
(Е11 + кА11)[(Е11 — кАи)У 1 — кАиУ 2 + кУ 1] + кАпй 2 = Ух + кУ 1,
или
(Е11 — к2 А2^ 1 + (Е11 + кА^—кА^ 2 + кУ 1] + кА1и 2 = V1 + кУ1,
или
(Е11 — к2 А121 — Еп)¥ 1 — (Е11 + кД 1)кА12 V 2 = —кА12 и 2 + к( Е11 — Е11 — кА11 )У 1
или
— к2 А2 V1 = (Е11 + кА^кА^ 2 — кА12и 2 — к2 А11У1,
или
кА121 V1 = — (Ец + кАп)А12V2 + А12и2 + кАцУ 1, (3.1.13) или
Г1 = — к ч( А121)—1[(Е11 + кА11) А12Г2 + А12и 2 + кА11У1]. (4.1.13)
Теперь подставим полученное выражение (4.1.13) в (4.1.11.2).
и 2 = А21 (А2)—1[( Е11 + кАп) А^ 2 + А12и 2 + кАцУ 1] + (Е22 — кА^У 2 + кУ 2, или
U2 = [ A21( A121) 1(E11 + hA11) A12 + (E22 - hA22)]V2 + A21( A121) 1 A12U2 +
+ h[ A21( A121) 1 A11Y1 + Y2 ] ,
или
[E22 - A21( A121) 1 A12]U2 = [ A21( A121) 1(E11 + hA11) A12 + (E22 - hA22)]V2 +
+ h[ A21( A121) 1 A11Y1 + Y2 ] .
(4.1.14)
Выводы
1). Получена новая вычислительная схема, в которой вычисляются только m << n коэффициентов вектора X(t) .
2). После перемножения всех матриц в (3.1.14) итоговое выражение будет иметь
вид
Bi( X 2)1 +1 = [B2 h + B3]( X 2)1 + h[ B4(Y i)i + (Y 2)1 ], (4.1.15)
где B1, B2, B3 - числовые (m x m) - матрицы, B4 - числовая (m x n) -
матрица. В практических задачах вектор Y(t) содержит в себе всего лишь несколько ненулевых коэффициентов. Поэтому на практике выражение (4.1.15) будет иметь более простой вид
Bi( X 2)1 +1 = [Bh + B3KX2X + h(Y 2)1 ]. (4.1.16)
3). По построению матрица B1 не вырождена и не зависит от параметра h.
Поэтому от нее может выть вычислена обратная матрица, после умножения на которую выражения (4.1.16), получим вычислительную схему в виде
(X2)i+1 = [Cih + C2KX2X + hC,{Y2\]. (4.1.17)
4). Трудоемкость вычислений по схеме (4.1.17) при измененном шаге составит
2 2 3m ВМО, а при постоянном шаге - m ВМО. При n = 500 и m = 10 для плотных
матриц скорость вычислений по сравнению с неявным методом Эйлера увеличится
примерно в 417000 раз, а по сравнению с явным методом Эйлера в 833 раза.
Список литературы
1. Новиков Л.С., Бабкин Г.В., Морозов Е.П., Колосов С. А., Крупников К.К., Милеев В.Н., Саенко В.С. Комплексная методология определения параметров электростатической зарядки, электрических полей и пробоев на космических аппаратах в условиях радиационной электризации. Руководство для конструкторов. - ЦНИИмаш, 1995. - 160 с.
2. Бабкин Г.В., Морозов Е.П. Активная защита космических аппаратов от статического электричества в орбитальных условиях: Справочное руководство для конструкторов. Королев. М.О.: ЦНИИмаш, 2000. 285 с.
3. Соколов А.Б., Марченков К.В., Борисов Н.И., Саенко В.С., Бабкин Г.В. Разработка алгоритма формирования структурной электрофизической модели космического аппарата на основе электрических схем, состоящих из фазовых параметрических макромоделей. // Космонавтика и ракетостроение. - Москва, ЦНИИмаш. - Вып. 3(52) - 2008. - с. 161-174.
4. Востриков А.В. Приближенный метод расчета растекания токов по элементам конструкции космического аппарата при электростатических разрядах. Технологии ЭМС №2(33). Москва 2010, с. 75-79.
5. Чуа Л.О., Лин Пен-Мин Машинный анализ электронных схем: Алгоритмы и вычислительные методы. Пер. с англ. - М.: Энергия, 1980. - 640 с.
