РАЗРАБОТКА МОДЕЛИ РАСПРЕДЕЛЕНИЯ ПЛОТНОСТИ ТОКОВ ПРИ ВОЗБУЖДЕНИИ ИОНОСФЕРЫ ВЫСОКОЧАСТОТНЫМ ОБЛУЧЕНИЕМ
ЭЛЕКТРОННЫЕ И ЭЛЕКТРОМАГНИТНЫЕ УСТРОЙСТВА
УДК 51.74: 533.9: 537.5: 537.876
A. Ю. ГРИШЕНЦЕВ, А. Г. КОРОБЕЙНИКОВ
РАЗРАБОТКА МОДЕЛИ РАСПРЕДЕЛЕНИЯ ПЛОТНОСТИ ТОКОВ ПРИ ВОЗБУЖДЕНИИ ИОНОСФЕРЫ ВЫСОКОЧАСТОТНЫМ ОБЛУЧЕНИЕМ
Разработана модель объемного распределения плотности токов в ионосфере в период облучения высокочастотным радиоизлучением при возникновении переходных процессов и в устоявшихся режимах для решения прямой и обратной задачи. Рассмотрены методы численного решения модели и расчета ее параметров.
Ключевые слова: ионосфера, распространение электромагнитных колебаний, математическое моделирование.
Введение. Моделирование протекания в ионосфере динамических процессов, связанных с возбуждением электромагнитного поля, актуально для радиосвязи, метеорологии и ряда других научно-практических направлений [1—8]. Моделирование распространения радиоволн в ионосфере является сложной многопараметрической задачей, пока не имеющей универсального решения [1—11]. Задача усложняется особенностями плазменного состояния газа, пространственной и динамической неоднородностью ионосферы, неоднородностью геомагнитного поля, наличием собственных токов в ионосфере и др. Модель, рассмотренная в настоящей работе, предназначена для численного анализа и изучения переходных процессов, а также пространственного распределения плотности тока в ионосфере при искусственном облучении высокочастотным радиоизлучением.
Геометрическая структура модели. Предлагаемая модель построена путем разбиения некоторого выделенного объема атмосферы (включая тропосферу, стратосферу и ионосферу) на дискретные части цилиндрическо-кольцевой формы с дальнейшим замещением их эквивалентными электрическими RLC-ветвями (рис. 1), впоследствии объединенными в совокупную RLC-схему. Использование слоистой структуры модели позволяет моделировать движение электронного газа.
Рассмотрим один из возможных способов разбиения конусного участка атмосферы. Разделив участок на дискреты, получим коаксиальную многослойную структуру, внутри которой можно выделить объемы, обладающие относительной однородностью свойств (рис. 1).
Область ионосферы разделена по горизонтали на K коаксиальных слоев (от нуля до K–1) поверхностями сечения, имеющими форму усеченного сферой конуса с шаровым сегментом в основании и относительно малым углом при вершине (такие усеченные конусы, при некоторых допущениях, можно считать цилиндрами); по вертикали — на M слоев (от нуля до M–1) поверхностями сечения, имеющими форму шарового сегмента. Таким образом, ионосфера дробится на множество коаксиальных цилиндрических трубок общим числом NR=(M+1)K,
ИЗВ. ВУЗОВ. ПРИБОРОСТРОЕНИЕ. 2010. Т. 53, № 12
42 A. Ю. Гришенцев, А. Г. Коробейников
каждая из которых характеризуется радиусом rmk , толщиной ∆rmk и высотой hm, где
m=0, 1, 2, …, M–1 и k=0, 1, 2, ..., K–1. Область тропосферы и стратосферы возможно объеди-
нить в единый слой M и разделить на K–1 коаксиальных трубок высотой hM, которые сверху примыкают к соответствующим трубкам из слоя ионосферы M–1, а снизу опираются на по-
верхность Земли.
∆r00
r00
∆r01
r01
h01
01 0
K–1
1
Ионосфера
M–1
hm c
c E
Земля
Рис. 1
Выбор значений hm для последующего моделирования принимается с учетом изменения следующих величин как функций высоты: значение электронной плотности Ne, эффективного числа упругих столкновений νeff электронов с нейтральными частицами и ионами, температуры (в том числе температурных инверсий), химического состава и плотности. Значения ∆rmk и Rmk выбираются так, чтобы можно было оценить неравномерность распределения плотности тока по ионосфере, вызванную проявлением поверхностного эффекта, т.е. с учетом глубины
проникновения δ электромагнитной волны:
δ
=
λ 2π
=
2, ωµγ
(1)
где ω — круговая частота, µ — абсолютная магнитная проницаемость вещества, γ — удельная объемная электропроводность вещества, λ — длина волны. Наиболее значительные изменения с высотой (h) аргументов функции (1) будут наблюдаться для γ в зависимости от значений Ne и νeff [12]. При выборе ∆rmk требуется оценивать уточненное значение δ с наибольшими значениями γ для слоев ионосферы, в которых необходимо учитывать неравномерность распределения плотности тока. Окончательное уточнение геометрических параметров модели производится исходя из требуемой точности и доступных вычислительных ресурсов, а также с учетом особенностей расчета эквивалентных электрических параметров.
Эквивалентная схема замещения. Эквивалентная схема (рис. 2) получается в результате замены каждого дискретного элемента модели набором RLC с определенными параметрами, значения которых могут динамически меняться во времени. В данной схеме элементы имеют коаксиальную структуру, каждый (см. рис. 1) обладает собственной индуктивностью L(mK+k) и активным сопротивлением R(mK+k), где m=0, …, M–1, k=0, ..., K–1. Кроме того, присутствуют межслойные электрические взаимодействия благодаря емкостным и индуктивным связям, индуктивно в модели взаимосвязаны все катушки, между соседними слоями присутствуют емкостные связи C(m(K–1)+k). Помимо емкостных связей, обусловленных токами смещения, присутствуют горизонтальные межслойные связи, обусловленные токами проводимости R(m(K–1)+k). Такой подход позволяет рассматривать модель в виде схемы, состоя-
ИЗВ. ВУЗОВ. ПРИБОРОСТРОЕНИЕ. 2010. Т. 53, № 12
Разработка модели распределения плотности токов при возбуждении ионосферы облучением 43 щей из обычных RLC-элементов, и анализировать переходные процессы и устоявшиеся режимы стандартными методами электротехники.
Рис. 2
Индексы элементов (см. рис. 2) соответствуют индексам геометрической структуры модели. Сквозная нумерация элементов RLC позволяет получить схему, имеющую стандартные обозначения, что может существенно облегчить задачу численного моделирования, например, при использовании пакетов схемотехнического моделирования.
Для составления системы дифференциальных уравнений первого порядка каждый элемент цепи выделим в отдельную ветвь. При отнесении ветвей графа схемы к дереву или связям будем придерживаться правила получения нормального дерева и нормального подграфа связей. К ветвям дерева последовательно отнесем ветви с источниками ЭДС, затем — ветви с конденсаторами, ветви с резисторами и в последнюю очередь — ветви с катушками индуктивности. В качестве связей (хорд) сначала выделим источники тока (в приведенном варианте схемы отсутствуют), затем — индуктивные элементы, резистивные ветви и в последнюю очередь — ветви с конденсаторами. Граф, полученный таким образом, называют нормальным [13].
Общее число узлов в графе схемы [(2M + 1)K + 2] , общее число ветвей [K (4M +1) +1 − 2M ].
Рассмотрим алгоритм получения матриц, необходимых для машинного решения системы дифференциальных уравнений, описывающих токи и напряжения в ветвях схемы. Занесем последовательно все ветви в общую таблицу DATA в соответствии с правилами получения нормального дерева и подграфа связей. Таблица DATA должна содержать следующий
ИЗВ. ВУЗОВ. ПРИБОРОСТРОЕНИЕ. 2010. Т. 53, № 12
44 A. Ю. Гришенцев, А. Г. Коробейников
минимальный набор полей: порядковый номер строки (далее соответствует номеру ветви), номера узлов стока и истока, тип и номер элемента в соответствии с электрической схемой.
Далее путем перебора строк DATA получим матрицу A, строки которой соответствуют номерам узлов, а столбцы — номерам ветвей. Номера ветвей соответствуют номеру элемента в массиве DATA.
В каждом столбце полной матрицы A должно быть две единицы, одна (+1) соответствует начальному узлу ветви (истоку), вторая (–1) — конечному (стоку). Для пассивных элементов цепи не имеет принципиального значения положение стока (истока). Полученная матрица содержит одну избыточную строку, являющуюся линейной комбинацией других, поэтому одну строку необходимо исключить.
Размерность матрицы A [(2M +1)K + 2 −1]×[K (4M +1) +1− 2M ], где [(2M +1)K + 2 −1] —
число узлов, [K (4M +1) +1 − 2M ] — число ветвей.
Следующей операцией является деление полученной матрицы на матрицу дерева Aд и
матрицу связей Aс путем последовательного перебора в направлении, соответствующем на-
правлению заполнения матрицы A из таблицы DATA. При формировании матрицы Aд достаточно, чтобы хотя бы один из узлов нового элемента не содержался в списке узлов, перечис-
ленных ранее (если это условие не выполняется, то элемент переносится в матрицу Aс). В результате необходимо получить матрицу следующего вида:
A = Aд Aс .
(2)
С учетом того что таблица DATA была получена в соответствии с правилами формиро-
вания нормального графа, матрица Aд будет содержать набор данных, соответствующий нормальному графу. В полученных матрицах дерева Aд и связей Aс номера ветвей могут быть расположены не по порядку, что может вызывать некоторые неудобства при хранении в памяти ЭВМ и вычислениях, привести нумерацию к упорядоченному виду можно с помощью упорядоченного перебора и присвоения новых номеров ветвям. В таблице DATA для возможности взаимообратных преобразований необходимо добавить поле, содержащее новую нумерацию ветвей. В результате стандартных преобразований матриц получим систему дифференциальных уравнений, которая в совокупности с начальными условиями представляет задачу Коши, такую систему можно решать различными методами численного интегрирования. Точность и устойчивость решения будут во многом зависеть от выбранного шага интегрирования, при этом слишком малые значения могут существенно увеличить время расчетов. При формировании начальных условий необходимо учитывать следующие правила:
iL (t − 0) = iL (t + 0) и uC (t − 0) = uC (t + 0) , где iL — ток через индуктивности в зависимости от
времени t, uC — напряжение на конденсаторе в зависимости от времени t.
Расчет эквивалентных электрических параметров геометрической модели. Расчет эквивалентных значений активных сопротивлений при движении тока в вертикальном направлении (см. рис. 1) для дискретных элементов модели с номерами k и m можно осуществить по формуле:
RmK +k
=
1 γ km
hm 2πrmk ∆rmk
.
(3)
Межслойное сопротивление при движении тока в горизонтальном направлении можно
рассчитать следующим образом. Рассчитаем элементарное приращение межслойного актив-
ного сопротивления для одного (внутреннего) полуслоя:
dRm′ (K −1)+k +NR
1 =
γ mk
dr
2π(rmk + dr ) hm
,
далее запишем в интегральном виде:
ИЗВ. ВУЗОВ. ПРИБОРОСТРОЕНИЕ. 2010. Т. 53, № 12
Разработка модели распределения плотности токов при возбуждении ионосферы облучением 45
rmk
+
∆rmk 2
∫ dr
Rm′ (K −1)+k +NR
=
1 γ mk
rmk
⎛ ⎜
rmk
+
∆rmk 2
⎞
=
1 γ mk
∆rmk
,
2πhm (2rmk + ∆rmk )
⎟
∫2πhm ⎜ rmk +
dr ⎟
⎜⎜⎝ rmk ⎠⎟⎟
для внешнего слоя на основании аналогичных преобразований получаем
( )Rm′′ (K −1)+k+NR
=
1 γ m( k +1)
2πhm
∆rm( k +1) 2rm(k+1) + ∆rm(k +1)
.
Полное межслойное сопротивление можно вычислить как
Rm(K −1)+k +NR
1 =
γ mk
2πhm
∆rmk (2rmk +
∆rmk
)
+
1 γ m( k +1)
∆rm( k +1)
.
2πhm (2rm(k+1) + ∆rm(k+1) )
(4)
При решении задач на расчет индуктивности [14] принято рассматривать случай посто-
янного тока и низкой частоты, когда поперечные линейные размеры проводника меньше, или значительно меньше, глубины проникновения тока δ (1); случай весьма высокой частоты, ко-
гда поперечные линейные размеры проводника больше или равны δ. При машинном расчете значения индуктивности (собственные и взаимные) целесообразно хранить в двумерном,
квадратном, массиве IND размером [(M +1)K ]× [(M + 1)K ]; элементы с одинаковыми номе-
рами строк и столбцов будут образовывать множество собственных индуктивностей, осталь-
ные элементы будут соответствовать взаимным индуктивностям, причем IND[i, j] = IND[ j,i] .
Рассмотрим некоторые примеры, используемые в модели для расчета индуктивностей.
В общем случае собственную индуктивность коаксиального проводника приближенно можно определить как
IND[mK + k, mK
+ k ] = LmK +k
=
µ0 hm 2π
(ln
2hm
rmk
+
∆rmk 2
−1− µ µ0
ln c) ,
(5)
где µ0 — магнитная постоянная, с — табличная константа [14], зависящая от отношения
внутреннего радиуса
r−
∆r 2
к внешнему
r+
∆r 2
.
В пределах одного горизонтального слоя взаимная индуктивность рассчитывается для
коаксиальных проводников в случае низкой частоты:
( ) ( )⎧ ⎡
⎛
⎪⎪M ⎨ ⎪
=
µ0 2π
⎢⎢ln ⎢ ⎣
q p
+
µ µ0
⎜ ⎜ ⎜⎝
s4 s2 −q2
ln
s q
+
n4 p2 − n2
2
ln
p n
−
1 2
⎝⎜⎜⎛1+
q2 s2 −q2
+
n2 p2 − n2
⎞ ⎠⎟⎟
⎞⎤ ⎟⎥ ⎟⎥ ⎟⎠⎥⎦
,
⎪⎩IND[i, j] = IND[ j, i] = M ,
(6)
где
n
=
ri
−
∆ri 2
,
p
=
ri
+
∆ri 2
,
q
=
rj
−
∆r j 2
,
s = rj
+
∆rj 2
, причем
rj
> ri , а
p = q.
Расчет взаимной индуктивности между M элементами из различных горизонтальных
слоев как между полыми цилиндрическими проводниками с общей осевой линией рассмот-
рен в [14], в общем случае значение M является функцией от линейных размеров
ri , ∆ri , hi , rj , ∆rj , hj дискретных проводников i, j и расстоянием между их центрами x:
ИЗВ. ВУЗОВ. ПРИБОРОСТРОЕНИЕ. 2010. Т. 53, № 12
46 A. Ю. Гришенцев, А. Г. Коробейников
⎪⎧M = f (ri , ∆ri , hi , rj , ∆rj , hj , x),
⎩⎨⎪IND[i, j] = IND[ j,i] = M .
(7)
Рассчитаем межслойную электрическую емкость Cmk как совокупность двух последо-
вательно включенных емкостей
Cm′ k
= 2πε0εmk
hm
ln
⎛ ⎜ ⎝
rmk
+ ∆rmk rmk
⎞ ⎟ ⎠
и
Cm′′k
=
2πε0εmk+1
hm
ln
⎜⎜⎜⎛⎜⎝
rmk
+1 + ∆rmk rmk +1
+1
⎠⎞⎟⎟⎟⎟
,
Cm′ k — емкость внутреннего полуслоя, Cm′′k — емкость внешнего полуслоя по отношению к
площади, разделяющей слои:
Cmk
=
Cm′ k Cm′′ k Cm′ k + Cm′′k
.
(8)
Предложенные способы расчета эквивалентных электрических параметров модели мо-
гут варьироваться в зависимости от требований к точности и особенностей численной реали-
зации модели.
Заключение. Разработанная модель применима к случаю искусственного возбуждения
ионосферы при помощи узконаправленного потока высокочастотного радиоизлучения. В ре-
зультате моделирования возможно получить трехмерную картину пространственного распре-
деления токов в ионосфере, определить динамику переходных процессов при изменении ин-
тенсивности излучения. Предложенная модель позволяет изучать процессы переизлучения и
модуляции радиоволн ионосферой, т.е. обратную задачу. На основе предложенной модели
возможно изучать некоторые особенности управления модуляцией сигналов при принуди-
тельном радиооблучении ионосферы.
СПИСОК ЛИТЕРАТУРЫ
1. Rycroft M. J. Solar-terrestrial energy program: handbook of ionospheric models // J. Atmospheric and SolarTerrestrial Physics. Elsevier Science Publishing Company, Inc., 1998. P. 403—404.
2. Guide to Reference and Standard Ionosphere Models. American Institute of Aeronautics and Astronautics, 1999. 55 p.
3. Rasmussen C. E., Schunk R. W. A three-dimensional time-dependent model of the plasmasphere // J. Geophys. Res. 1990. Vol. 95. P. 6133—6144.
4. Смирнов В. М., Смирнова Е. В. Реконструкция пространственно-временной структуры ионосферы по данным спутниковых наблюдений // Современные проблемы дистанционного зондирования Земли из космоса. 2008. Т. 5. С. 561—566.
5. Думин Ю. В. Концентрация носителей заряда в метастабильной плазме с сильной кулоновской неидеальностью // Прикладная физика. 1999. № 5. С. 18—21.
6. Richmond A. D., Ridley E. C., Roble R. G. A thermosphere/ionosphere general circulation model with coupled electrodynamics // Geophys. Res. Lett. 1992. Vol. 19, N 6. P. 601—604.
7. Голыгин В. А., Сажин В. И., Унучков В. Е. Коррекция модели ионосферы по данным о максимальноприменимых частотах реперных радиолиний // Исследовано в России. 2006 [Электронный ресурс]: .
8. Намгаладзе А. А., Юрик Р. Ю. Математическое моделирование возмущений верхней атмосферы Земли [Электронный ресурс]: .
9. Намгаладзе А. А., Мартыненко О. В., Волков М. А., Намгаладзе А. Н., Юрик Р. Ю. Математическое моделирование крупномасштабных возмущений верхней атмосферы Земли // Моделирование процессов в верхней полярной атмосфере. Апатиты: ПГИ КНЦ РАН, 1998. С. 167—249.
10. Гусятинский И. А., Немировский А. С., Соколов А. В., Троицкий В. Н. Дальняя тропосферная радиосвязь. М.: Связь, 1968. 248 с.
ИЗВ. ВУЗОВ. ПРИБОРОСТРОЕНИЕ. 2010. Т. 53, № 12
Дискретный интегродифференцирующий фильтр
47
11. Хлыбов Е. С., Гаврилов Б. Г., Егоров Д. Е. Модельные исследования ионосферных токов, электрических полей и ускорений заряженных частиц // Сб. тр. XLVII науч. конф. МФТИ. М.: Изд-во МФТИ, 2002.
12. Калинин А. И., Черенкова Е. Л. Распространение радиоволн и работа радиолиний. М.: Связь, 1971. 450 с.
13. Демерчан К. С., Нейман Л. Р., Коровкин Н. В., Чечурин В. Л. Теоретические основы электротехники. СПб: Питер, 2006. 576 с.
14. Калантаров П. Л., Цейтлин Л. А. Расчет индуктивностей. Л.: Энергоатомиздат, 1986. 488 c.
Алексей Юрьевич Гришенцев Анатолий Григорьевич Коробейников
Сведения об авторах — канд. техн. наук; Санкт-Петербургский государственный универси-
тет информационных технологий, механики и оптики, кафедра проектирования компьютерных систем; E-mail: tigerpost@ya.ru — д-р техн. наук, профессор; Санкт-Петербургский государственный университет информационных технологий, механики и оптики, кафедра проектирования компьютерных систем
Рекомендована кафедрой проектирования компьютерных систем
Поступила в редакцию 16.03.10 г.
ИЗВ. ВУЗОВ. ПРИБОРОСТРОЕНИЕ. 2010. Т. 53, № 12
УДК 51.74: 533.9: 537.5: 537.876
A. Ю. ГРИШЕНЦЕВ, А. Г. КОРОБЕЙНИКОВ
РАЗРАБОТКА МОДЕЛИ РАСПРЕДЕЛЕНИЯ ПЛОТНОСТИ ТОКОВ ПРИ ВОЗБУЖДЕНИИ ИОНОСФЕРЫ ВЫСОКОЧАСТОТНЫМ ОБЛУЧЕНИЕМ
Разработана модель объемного распределения плотности токов в ионосфере в период облучения высокочастотным радиоизлучением при возникновении переходных процессов и в устоявшихся режимах для решения прямой и обратной задачи. Рассмотрены методы численного решения модели и расчета ее параметров.
Ключевые слова: ионосфера, распространение электромагнитных колебаний, математическое моделирование.
Введение. Моделирование протекания в ионосфере динамических процессов, связанных с возбуждением электромагнитного поля, актуально для радиосвязи, метеорологии и ряда других научно-практических направлений [1—8]. Моделирование распространения радиоволн в ионосфере является сложной многопараметрической задачей, пока не имеющей универсального решения [1—11]. Задача усложняется особенностями плазменного состояния газа, пространственной и динамической неоднородностью ионосферы, неоднородностью геомагнитного поля, наличием собственных токов в ионосфере и др. Модель, рассмотренная в настоящей работе, предназначена для численного анализа и изучения переходных процессов, а также пространственного распределения плотности тока в ионосфере при искусственном облучении высокочастотным радиоизлучением.
Геометрическая структура модели. Предлагаемая модель построена путем разбиения некоторого выделенного объема атмосферы (включая тропосферу, стратосферу и ионосферу) на дискретные части цилиндрическо-кольцевой формы с дальнейшим замещением их эквивалентными электрическими RLC-ветвями (рис. 1), впоследствии объединенными в совокупную RLC-схему. Использование слоистой структуры модели позволяет моделировать движение электронного газа.
Рассмотрим один из возможных способов разбиения конусного участка атмосферы. Разделив участок на дискреты, получим коаксиальную многослойную структуру, внутри которой можно выделить объемы, обладающие относительной однородностью свойств (рис. 1).
Область ионосферы разделена по горизонтали на K коаксиальных слоев (от нуля до K–1) поверхностями сечения, имеющими форму усеченного сферой конуса с шаровым сегментом в основании и относительно малым углом при вершине (такие усеченные конусы, при некоторых допущениях, можно считать цилиндрами); по вертикали — на M слоев (от нуля до M–1) поверхностями сечения, имеющими форму шарового сегмента. Таким образом, ионосфера дробится на множество коаксиальных цилиндрических трубок общим числом NR=(M+1)K,
ИЗВ. ВУЗОВ. ПРИБОРОСТРОЕНИЕ. 2010. Т. 53, № 12
42 A. Ю. Гришенцев, А. Г. Коробейников
каждая из которых характеризуется радиусом rmk , толщиной ∆rmk и высотой hm, где
m=0, 1, 2, …, M–1 и k=0, 1, 2, ..., K–1. Область тропосферы и стратосферы возможно объеди-
нить в единый слой M и разделить на K–1 коаксиальных трубок высотой hM, которые сверху примыкают к соответствующим трубкам из слоя ионосферы M–1, а снизу опираются на по-
верхность Земли.
∆r00
r00
∆r01
r01
h01
01 0
K–1
1
Ионосфера
M–1
hm c
c E
Земля
Рис. 1
Выбор значений hm для последующего моделирования принимается с учетом изменения следующих величин как функций высоты: значение электронной плотности Ne, эффективного числа упругих столкновений νeff электронов с нейтральными частицами и ионами, температуры (в том числе температурных инверсий), химического состава и плотности. Значения ∆rmk и Rmk выбираются так, чтобы можно было оценить неравномерность распределения плотности тока по ионосфере, вызванную проявлением поверхностного эффекта, т.е. с учетом глубины
проникновения δ электромагнитной волны:
δ
=
λ 2π
=
2, ωµγ
(1)
где ω — круговая частота, µ — абсолютная магнитная проницаемость вещества, γ — удельная объемная электропроводность вещества, λ — длина волны. Наиболее значительные изменения с высотой (h) аргументов функции (1) будут наблюдаться для γ в зависимости от значений Ne и νeff [12]. При выборе ∆rmk требуется оценивать уточненное значение δ с наибольшими значениями γ для слоев ионосферы, в которых необходимо учитывать неравномерность распределения плотности тока. Окончательное уточнение геометрических параметров модели производится исходя из требуемой точности и доступных вычислительных ресурсов, а также с учетом особенностей расчета эквивалентных электрических параметров.
Эквивалентная схема замещения. Эквивалентная схема (рис. 2) получается в результате замены каждого дискретного элемента модели набором RLC с определенными параметрами, значения которых могут динамически меняться во времени. В данной схеме элементы имеют коаксиальную структуру, каждый (см. рис. 1) обладает собственной индуктивностью L(mK+k) и активным сопротивлением R(mK+k), где m=0, …, M–1, k=0, ..., K–1. Кроме того, присутствуют межслойные электрические взаимодействия благодаря емкостным и индуктивным связям, индуктивно в модели взаимосвязаны все катушки, между соседними слоями присутствуют емкостные связи C(m(K–1)+k). Помимо емкостных связей, обусловленных токами смещения, присутствуют горизонтальные межслойные связи, обусловленные токами проводимости R(m(K–1)+k). Такой подход позволяет рассматривать модель в виде схемы, состоя-
ИЗВ. ВУЗОВ. ПРИБОРОСТРОЕНИЕ. 2010. Т. 53, № 12
Разработка модели распределения плотности токов при возбуждении ионосферы облучением 43 щей из обычных RLC-элементов, и анализировать переходные процессы и устоявшиеся режимы стандартными методами электротехники.
Рис. 2
Индексы элементов (см. рис. 2) соответствуют индексам геометрической структуры модели. Сквозная нумерация элементов RLC позволяет получить схему, имеющую стандартные обозначения, что может существенно облегчить задачу численного моделирования, например, при использовании пакетов схемотехнического моделирования.
Для составления системы дифференциальных уравнений первого порядка каждый элемент цепи выделим в отдельную ветвь. При отнесении ветвей графа схемы к дереву или связям будем придерживаться правила получения нормального дерева и нормального подграфа связей. К ветвям дерева последовательно отнесем ветви с источниками ЭДС, затем — ветви с конденсаторами, ветви с резисторами и в последнюю очередь — ветви с катушками индуктивности. В качестве связей (хорд) сначала выделим источники тока (в приведенном варианте схемы отсутствуют), затем — индуктивные элементы, резистивные ветви и в последнюю очередь — ветви с конденсаторами. Граф, полученный таким образом, называют нормальным [13].
Общее число узлов в графе схемы [(2M + 1)K + 2] , общее число ветвей [K (4M +1) +1 − 2M ].
Рассмотрим алгоритм получения матриц, необходимых для машинного решения системы дифференциальных уравнений, описывающих токи и напряжения в ветвях схемы. Занесем последовательно все ветви в общую таблицу DATA в соответствии с правилами получения нормального дерева и подграфа связей. Таблица DATA должна содержать следующий
ИЗВ. ВУЗОВ. ПРИБОРОСТРОЕНИЕ. 2010. Т. 53, № 12
44 A. Ю. Гришенцев, А. Г. Коробейников
минимальный набор полей: порядковый номер строки (далее соответствует номеру ветви), номера узлов стока и истока, тип и номер элемента в соответствии с электрической схемой.
Далее путем перебора строк DATA получим матрицу A, строки которой соответствуют номерам узлов, а столбцы — номерам ветвей. Номера ветвей соответствуют номеру элемента в массиве DATA.
В каждом столбце полной матрицы A должно быть две единицы, одна (+1) соответствует начальному узлу ветви (истоку), вторая (–1) — конечному (стоку). Для пассивных элементов цепи не имеет принципиального значения положение стока (истока). Полученная матрица содержит одну избыточную строку, являющуюся линейной комбинацией других, поэтому одну строку необходимо исключить.
Размерность матрицы A [(2M +1)K + 2 −1]×[K (4M +1) +1− 2M ], где [(2M +1)K + 2 −1] —
число узлов, [K (4M +1) +1 − 2M ] — число ветвей.
Следующей операцией является деление полученной матрицы на матрицу дерева Aд и
матрицу связей Aс путем последовательного перебора в направлении, соответствующем на-
правлению заполнения матрицы A из таблицы DATA. При формировании матрицы Aд достаточно, чтобы хотя бы один из узлов нового элемента не содержался в списке узлов, перечис-
ленных ранее (если это условие не выполняется, то элемент переносится в матрицу Aс). В результате необходимо получить матрицу следующего вида:
A = Aд Aс .
(2)
С учетом того что таблица DATA была получена в соответствии с правилами формиро-
вания нормального графа, матрица Aд будет содержать набор данных, соответствующий нормальному графу. В полученных матрицах дерева Aд и связей Aс номера ветвей могут быть расположены не по порядку, что может вызывать некоторые неудобства при хранении в памяти ЭВМ и вычислениях, привести нумерацию к упорядоченному виду можно с помощью упорядоченного перебора и присвоения новых номеров ветвям. В таблице DATA для возможности взаимообратных преобразований необходимо добавить поле, содержащее новую нумерацию ветвей. В результате стандартных преобразований матриц получим систему дифференциальных уравнений, которая в совокупности с начальными условиями представляет задачу Коши, такую систему можно решать различными методами численного интегрирования. Точность и устойчивость решения будут во многом зависеть от выбранного шага интегрирования, при этом слишком малые значения могут существенно увеличить время расчетов. При формировании начальных условий необходимо учитывать следующие правила:
iL (t − 0) = iL (t + 0) и uC (t − 0) = uC (t + 0) , где iL — ток через индуктивности в зависимости от
времени t, uC — напряжение на конденсаторе в зависимости от времени t.
Расчет эквивалентных электрических параметров геометрической модели. Расчет эквивалентных значений активных сопротивлений при движении тока в вертикальном направлении (см. рис. 1) для дискретных элементов модели с номерами k и m можно осуществить по формуле:
RmK +k
=
1 γ km
hm 2πrmk ∆rmk
.
(3)
Межслойное сопротивление при движении тока в горизонтальном направлении можно
рассчитать следующим образом. Рассчитаем элементарное приращение межслойного актив-
ного сопротивления для одного (внутреннего) полуслоя:
dRm′ (K −1)+k +NR
1 =
γ mk
dr
2π(rmk + dr ) hm
,
далее запишем в интегральном виде:
ИЗВ. ВУЗОВ. ПРИБОРОСТРОЕНИЕ. 2010. Т. 53, № 12
Разработка модели распределения плотности токов при возбуждении ионосферы облучением 45
rmk
+
∆rmk 2
∫ dr
Rm′ (K −1)+k +NR
=
1 γ mk
rmk
⎛ ⎜
rmk
+
∆rmk 2
⎞
=
1 γ mk
∆rmk
,
2πhm (2rmk + ∆rmk )
⎟
∫2πhm ⎜ rmk +
dr ⎟
⎜⎜⎝ rmk ⎠⎟⎟
для внешнего слоя на основании аналогичных преобразований получаем
( )Rm′′ (K −1)+k+NR
=
1 γ m( k +1)
2πhm
∆rm( k +1) 2rm(k+1) + ∆rm(k +1)
.
Полное межслойное сопротивление можно вычислить как
Rm(K −1)+k +NR
1 =
γ mk
2πhm
∆rmk (2rmk +
∆rmk
)
+
1 γ m( k +1)
∆rm( k +1)
.
2πhm (2rm(k+1) + ∆rm(k+1) )
(4)
При решении задач на расчет индуктивности [14] принято рассматривать случай посто-
янного тока и низкой частоты, когда поперечные линейные размеры проводника меньше, или значительно меньше, глубины проникновения тока δ (1); случай весьма высокой частоты, ко-
гда поперечные линейные размеры проводника больше или равны δ. При машинном расчете значения индуктивности (собственные и взаимные) целесообразно хранить в двумерном,
квадратном, массиве IND размером [(M +1)K ]× [(M + 1)K ]; элементы с одинаковыми номе-
рами строк и столбцов будут образовывать множество собственных индуктивностей, осталь-
ные элементы будут соответствовать взаимным индуктивностям, причем IND[i, j] = IND[ j,i] .
Рассмотрим некоторые примеры, используемые в модели для расчета индуктивностей.
В общем случае собственную индуктивность коаксиального проводника приближенно можно определить как
IND[mK + k, mK
+ k ] = LmK +k
=
µ0 hm 2π
(ln
2hm
rmk
+
∆rmk 2
−1− µ µ0
ln c) ,
(5)
где µ0 — магнитная постоянная, с — табличная константа [14], зависящая от отношения
внутреннего радиуса
r−
∆r 2
к внешнему
r+
∆r 2
.
В пределах одного горизонтального слоя взаимная индуктивность рассчитывается для
коаксиальных проводников в случае низкой частоты:
( ) ( )⎧ ⎡
⎛
⎪⎪M ⎨ ⎪
=
µ0 2π
⎢⎢ln ⎢ ⎣
q p
+
µ µ0
⎜ ⎜ ⎜⎝
s4 s2 −q2
ln
s q
+
n4 p2 − n2
2
ln
p n
−
1 2
⎝⎜⎜⎛1+
q2 s2 −q2
+
n2 p2 − n2
⎞ ⎠⎟⎟
⎞⎤ ⎟⎥ ⎟⎥ ⎟⎠⎥⎦
,
⎪⎩IND[i, j] = IND[ j, i] = M ,
(6)
где
n
=
ri
−
∆ri 2
,
p
=
ri
+
∆ri 2
,
q
=
rj
−
∆r j 2
,
s = rj
+
∆rj 2
, причем
rj
> ri , а
p = q.
Расчет взаимной индуктивности между M элементами из различных горизонтальных
слоев как между полыми цилиндрическими проводниками с общей осевой линией рассмот-
рен в [14], в общем случае значение M является функцией от линейных размеров
ri , ∆ri , hi , rj , ∆rj , hj дискретных проводников i, j и расстоянием между их центрами x:
ИЗВ. ВУЗОВ. ПРИБОРОСТРОЕНИЕ. 2010. Т. 53, № 12
46 A. Ю. Гришенцев, А. Г. Коробейников
⎪⎧M = f (ri , ∆ri , hi , rj , ∆rj , hj , x),
⎩⎨⎪IND[i, j] = IND[ j,i] = M .
(7)
Рассчитаем межслойную электрическую емкость Cmk как совокупность двух последо-
вательно включенных емкостей
Cm′ k
= 2πε0εmk
hm
ln
⎛ ⎜ ⎝
rmk
+ ∆rmk rmk
⎞ ⎟ ⎠
и
Cm′′k
=
2πε0εmk+1
hm
ln
⎜⎜⎜⎛⎜⎝
rmk
+1 + ∆rmk rmk +1
+1
⎠⎞⎟⎟⎟⎟
,
Cm′ k — емкость внутреннего полуслоя, Cm′′k — емкость внешнего полуслоя по отношению к
площади, разделяющей слои:
Cmk
=
Cm′ k Cm′′ k Cm′ k + Cm′′k
.
(8)
Предложенные способы расчета эквивалентных электрических параметров модели мо-
гут варьироваться в зависимости от требований к точности и особенностей численной реали-
зации модели.
Заключение. Разработанная модель применима к случаю искусственного возбуждения
ионосферы при помощи узконаправленного потока высокочастотного радиоизлучения. В ре-
зультате моделирования возможно получить трехмерную картину пространственного распре-
деления токов в ионосфере, определить динамику переходных процессов при изменении ин-
тенсивности излучения. Предложенная модель позволяет изучать процессы переизлучения и
модуляции радиоволн ионосферой, т.е. обратную задачу. На основе предложенной модели
возможно изучать некоторые особенности управления модуляцией сигналов при принуди-
тельном радиооблучении ионосферы.
СПИСОК ЛИТЕРАТУРЫ
1. Rycroft M. J. Solar-terrestrial energy program: handbook of ionospheric models // J. Atmospheric and SolarTerrestrial Physics. Elsevier Science Publishing Company, Inc., 1998. P. 403—404.
2. Guide to Reference and Standard Ionosphere Models. American Institute of Aeronautics and Astronautics, 1999. 55 p.
3. Rasmussen C. E., Schunk R. W. A three-dimensional time-dependent model of the plasmasphere // J. Geophys. Res. 1990. Vol. 95. P. 6133—6144.
4. Смирнов В. М., Смирнова Е. В. Реконструкция пространственно-временной структуры ионосферы по данным спутниковых наблюдений // Современные проблемы дистанционного зондирования Земли из космоса. 2008. Т. 5. С. 561—566.
5. Думин Ю. В. Концентрация носителей заряда в метастабильной плазме с сильной кулоновской неидеальностью // Прикладная физика. 1999. № 5. С. 18—21.
6. Richmond A. D., Ridley E. C., Roble R. G. A thermosphere/ionosphere general circulation model with coupled electrodynamics // Geophys. Res. Lett. 1992. Vol. 19, N 6. P. 601—604.
7. Голыгин В. А., Сажин В. И., Унучков В. Е. Коррекция модели ионосферы по данным о максимальноприменимых частотах реперных радиолиний // Исследовано в России. 2006 [Электронный ресурс]: .
8. Намгаладзе А. А., Юрик Р. Ю. Математическое моделирование возмущений верхней атмосферы Земли [Электронный ресурс]: .
9. Намгаладзе А. А., Мартыненко О. В., Волков М. А., Намгаладзе А. Н., Юрик Р. Ю. Математическое моделирование крупномасштабных возмущений верхней атмосферы Земли // Моделирование процессов в верхней полярной атмосфере. Апатиты: ПГИ КНЦ РАН, 1998. С. 167—249.
10. Гусятинский И. А., Немировский А. С., Соколов А. В., Троицкий В. Н. Дальняя тропосферная радиосвязь. М.: Связь, 1968. 248 с.
ИЗВ. ВУЗОВ. ПРИБОРОСТРОЕНИЕ. 2010. Т. 53, № 12
Дискретный интегродифференцирующий фильтр
47
11. Хлыбов Е. С., Гаврилов Б. Г., Егоров Д. Е. Модельные исследования ионосферных токов, электрических полей и ускорений заряженных частиц // Сб. тр. XLVII науч. конф. МФТИ. М.: Изд-во МФТИ, 2002.
12. Калинин А. И., Черенкова Е. Л. Распространение радиоволн и работа радиолиний. М.: Связь, 1971. 450 с.
13. Демерчан К. С., Нейман Л. Р., Коровкин Н. В., Чечурин В. Л. Теоретические основы электротехники. СПб: Питер, 2006. 576 с.
14. Калантаров П. Л., Цейтлин Л. А. Расчет индуктивностей. Л.: Энергоатомиздат, 1986. 488 c.
Алексей Юрьевич Гришенцев Анатолий Григорьевич Коробейников
Сведения об авторах — канд. техн. наук; Санкт-Петербургский государственный универси-
тет информационных технологий, механики и оптики, кафедра проектирования компьютерных систем; E-mail: tigerpost@ya.ru — д-р техн. наук, профессор; Санкт-Петербургский государственный университет информационных технологий, механики и оптики, кафедра проектирования компьютерных систем
Рекомендована кафедрой проектирования компьютерных систем
Поступила в редакцию 16.03.10 г.
ИЗВ. ВУЗОВ. ПРИБОРОСТРОЕНИЕ. 2010. Т. 53, № 12