Определение характеристик ИНЧ-волн, наиболее сильно реагирующих на незначительные изменения электронной плотности ионосферы в области высоких широт
Автор: Ахметов О.И., Мингалев И.В., Мингалев О.В., Суворова З.В., Белаховский В.Б., Черняков С.М.
Журнал: Солнечно-земная физика @solnechno-zemnaya-fizika
Статья в выпуске: 4 т.5, 2019 года.
Бесплатный доступ
Методами численного эксперимента исследуется задача определения характеристик электромагнитных волн УНЧ-диапазона (0.3-3 кГц), регистрируемых на уровне приземного слоя и несущих максимальное количество информации о состоянии волновода Земля-ионосфера. Проанализировано влияние горизонтальной пространственной структуры электронной плотности волновода Земля-ионосфера на особенности распространения электромагнитных волн. Выявлены характеристики, позволяющие регистрировать их инструментальными методами наблюдений в условиях слабо возмущенной ионосферы. Профили концентрации, используемые в численных экспериментах, получены по данным радара частичных отражений Полярного геофизического института, расположенного на радиофизическом полигоне «Туманный» Мурманской области (69.0° N, 35.7° E), с использованием модели IRI2016 во время солнечной вспышки 15 марта и последующей магнитной бури 17 марта 2013 г. Используемая в работе модель распространения электромагнитных сигналов разработана авторами...
Унч-волны, численное моделирование, ионосфера
Короткий адрес: https://sciup.org/142222503
IDR: 142222503 | УДК: 537.877, | DOI: 10.12737/szf-54201911
Determination of ELF-wave characteristics most strongly reacting to minor changes of ionospheric electron density in a high-latitude region
In this paper, we use numerical experiment methods to address the problem of determining characteristics of ULF (0.3-3 kHz) electromagnetic waves recorded in the surface layer and providing the maximum amount of information about the Earth-ionosphere waveguide. We have analyzed the effect of the horizontal spatial structure of electron density of the Earth-ionosphere waveguide on propagation of electromagnetic waves. We have identified characteristics that allow us to record them by instrumental methods in conditions of weakly disturbed ionosphere. The density profiles used in numerical experiments have been obtained from data acquired by the Partial Reflection Radar at the Polar Geophysical Institute, located at the radiophysical observatory Tumanny in the Murmansk Region (69.0° N, 35.7° E), and by the IRI2016 model during the March 15, 2013 solar flare and the subsequent magnetic storm on March 17, 2013. The electromagnetic signal propagation model used in this work is the adaptation of gas-hydrodynamic methods to electrodynamic applications. The model is based on the scheme of upwind approximation of spatial derivatives (Godunov’s method with correction of streams). We also use splitting by spatial directions and physical processes. Signal field attenuation due to conductivity and its rotation due to Hall conductivity of the medium are considered in separate splitting steps by analytical formulas.
Текст научной статьи Определение характеристик ИНЧ-волн, наиболее сильно реагирующих на незначительные изменения электронной плотности ионосферы в области высоких широт
Одними из основных задач радиофизики в области распространения электромагнитных сигналов в атмосфере Земли являются поиск методов обеспечения стабильной связи или радиолокации на различных частотах; определение условий, когда эта связь невозможна; прогнозирование исключения таких условий и поиск мешающих связи или радиолокации факторов. Диагностика D-области ионосферы важна для многих физических приложений. В настоящее время основным отработанным средством такой диагностики являются радары частичных отражений. Ионозонды и радары некогерентного рассеяния формально захватывают D-область, но имеют там точность хуже, чем радары частичных отражений [Акимов, Калинин, 2017] . Состояние D-слоя ионосферы влияет на характеристики электромагнитных сигналов в диапазонах ИНЧ (0.3–3 кГц) и ОНЧ (3–30 кГц). Однако остается вопрос извлечения максимального количества информации о состоянии D-слоя из наземных измерений электромагнитных сигналов в указанном диапазоне.
В работе с помощью численного моделирования изучается влияние небольших (3–5 %) изменений концентрации электронов в D-слое ионосферы на характеристики электромагнитного сигнала на частоте 1500 Гц. Небольшие изменения концентрации электронов выбраны авторами с целью определения границ возможности диагностики D-слоя ионосферы сигналами ИНЧ-диапазона. Показано, что эти изменения не влияют на амплитуду отдельных компонент электромагнитного поля сигнала на уровне поверхности Земли, оказывают незначительное влияние на отношение амплитуд напряженностей электрического и магнитного полей и заметное влияние на разность фаз между напряженностями электрического и магнитного полей. Также исследовано влияние горизонтальных размеров области повышенной концентрации на перечисленные характеристики сигнала.
Для численного моделирования распространения сигналов ИНЧ-диапазона авторы использовали ранее разработанную ими модель распространения электромагнитных сигналов в системе литосфера— атмосфера—ионосфера. В модели численно интегрируется по времени система уравнений Максвелла на регулярной пространственной сетке. Для этого используется явная схема расщепления по пространственным направлениям и физическим процессам с противопотоковой аппроксимацией пространственных производных (метод Годунова с коррекцией потоков). Эта схема является консервативной, монотонной, имеет второй порядок точности по времени и третий по пространственным переменным.
В представленных численных экспериментах авторы используют прямоугольную область моделирования и декартову систему координат, пренебрегая кривизной сферических слоев литосферы, атмосферы и ионосферы по следующим причинам. Наибольший из горизонтальных размеров области моделирования примерно в шесть раз меньше радиуса Земли, а длина волны (~200 км) сравнима с размерами области моделирования. При таком соотношении масштабов сигнал распространяется горизонтально без многократных отражений от ионосферы и литосферы. Целью моделирования является исследование влияния локальных неоднородностей концентрации электронов в D-области ионосферы на характеристики электромагнитного сигнала. Плоская геометрия области моделирования для всех представленных численных экспериментов будет давать одинаковую малую ошибку, которая не окажет влияния на обсуждаемые закономерности в характеристиках сигнала.
ПОСТАНОВКА ЗАДАЧИ И ЗАДАНИЕ СРЕДЫ
В представленной работе авторы подробно рассматривают влияние незначительного повышения электронной плотности в нижней ионосфере в области высоких широт на прохождение электромагнитных волн ИНЧ-диапазона при различной конфигурации горизонтальных профилей электронной концентрации. Параметры среды задавались следующим образом: концентрация нейтральных частиц — с помощью эмпирической модели NRLMSISE-00; концентрация и температура электронов — с помощью эмпирической модели IRI 2016 с коррекцией по записям средневолнового радиолокатора вертикального излучения для исследования нижней ионосферы, расположенного на радиофизическом полигоне «Туманный» Полярного геофизического института в Мурманской области (69.0° N, 35.7° E) [Терещенко и др., 2003] , во время и после солнечной вспышки 15.03.2013, а именно 15.03.2013 в 09:00 UT (возмущенные условия, вспышка), 16.03.2013 в 09:00 UT (спокойные условия), 17.03.2013 в 06:00 UT (возмущенные условия, магнитная буря) и 16.03.2013 в 06:00 UT (спокойные условия). Вертикальные профили электронной плотности в указанных условиях показаны на рис. 1. Коррекция заключалась в перенесении тонкой структуры D- и E-областей, полученной по данным радиолокатора, на профили, полученные с помощью эмпирических моделей IRI 2016 и NRLMSISE.
Магнитная буря 17 марта 2013 г. была вызвана солнечной вспышкой М-класса 15 марта и началась с приходом к магнитосфере межпланетной ударной волны, которая проявилась на земной поверхности
Рис. 1. Аппроксимированные высотные профили концентрации электронов 15.03.2013 в 09:00 UT ( а ); 16.03.2013 в 09:00 UT ( б ); 17.03.2013 в 06:00 UT ( в ) и 16.03.2013 в 06:00 UT ( г )
Рис. 2. Горизонтальные профили концентрации электронов для эксперимента по распространению электромагнитной волны из спокойной области в возмущенную ( а ); из возмущенной области в спокойную ( б ); сквозь полосу повышенной концентрации электронов ( в )
в виде SSC (sudden storm commencement) в 06:00 UT. В этот момент скорость солнечного ветра резко увеличилась с ~400 до ~650–700 км/с, B z -компонента ММП стала отрицательной, что обеспечило постоянное поступление энергии в магнитосферу. Индекс SYM-H , характеризующий интенсивность магнитной бури, упал до –100 нТл и оставался на этом уровне. Авроральный AE -индекс резко возрос до ~1000 нТл и оставался повышенным. АЕ -индекс показывал в 16:00 UT появление другой авроральной активизации (увеличение AE до ~2500 нТл).
Несмотря на высокий класс вспышки, радар частичных отражений на полигоне «Туманный» и модель IRI показали незначительное повышение концентрации электронов на высотах D–E-областей ионосферы во время последующей магнитной бури. Авторы полагают, что такие условия хорошо подходят для разделения технически наблюдаемых и ненаблюдаемых характеристик электромагнитных волн и выделения тех из них, которые могут быть использованы для создания методов косвенных исследований ионосферы.
Профиль проводимости литосферы был взят из работы [Korja et al., 2002] , проводимость у поверхности составляла 2∙10–5 См.
Пространственная горизонтальная структура электронной плотности тоже может оказывать влияние на получаемые результаты. Для обнаружения этого влияния были проведены три серии экспериментов с горизонтальными профилями электронной плотности разных типов. Первая серия моделирует переход из спокойной области в область повышенной концентрации электронов. Горизонтальный профиль электронной концентрации показан на рис. 2, а. Вторая серия моделирует переход из области повышенной концентрации в спокойную область (профиль на рис. 2, б). Третья серия — прохождение полосы повышенной концентрации (рис. 2, в). Для всех серий переход осуществлялся по функции Гаусса. При этом в первых двух изменялся горизонтальный градиент для обнаружения его влияния в области перехода на распространение электромагнитных волн, в третьей изменялась ширина полосы с 0.25λ до λ, где λ — длина волны. Магнитное поле было взято вертикальным, что характерно для высоких широт, его индукция составляла 5.3^10-5 нТл.
ОБЛАСТЬ МОДЕЛИРОВАНИЯ, ПРОСТРАНСТВЕННАЯ СЕТКА И ИСТОЧНИК СИГНАЛА
Во всех поставленных в работе численных экспериментах область моделирования представляла собой параллелепипед с основанием 512x1280 км, высотой в атмосфере 150 км и глубиной в литосфере 50 км. Шаги сетки по горизонтали составляли 2 км, по вертикали в атмосфере 2 км и 0.5 км в литосфере. Шаг интегрирования по времени составлял 4^10-6 с.
На всех границах области моделирования применялись специальные условия свободного ухода волны для минимизации влияния нефизических отражений на результаты моделирования. Авторы применяли для моделирования постоянно действующий гармонический источник сигнала, а результаты анализировались при прекращении всех процессов, связанных с прохождением первичного фронта сигнала (35-й период колебаний). Поэтому, несмотря на почти 100 % отражение сигнала от высокопроводящих литосферы и ионосферы, находящихся внутри области моделирования, для исключения любой возможности накопления ошибок на вертикальных гранях расчетной области применялись такие же условия, как и на горизонтальных.
В качестве источника сигнала во всех представленных экспериментах используется поверхность одной из граней расчетной области. Такой источник сигнала позволяет задать не только амплитуду сигнала во времени, но и ее распределение в пространстве; задержками можно сформировать фронт волны необходимой формы подобно тому, как это делается в плоской эквидистантной фазированной антенной решетке [Воскресенский и др., 2012] . В экспериментах, обсуждаемых в данной работе, моделировалась плоская волна частотой 1500 Гц (длина волны ~200 км), излучаемая в область под прямым углом к плоскости источника.
ОПИСАНИЕ МОДЕЛИ
Используемая в представленной работе модель распространения электромагнитных сигналов в волноводе Земля—ионосфера с учетом инерции электронов в ионосфере построена на численном интегрировании по времени системы уравнений Максвелла на регулярной пространственной сетке с помощью явной схемы. Используется расщепление по пространственным направлениям и по физическим процессам. В схеме чередуются шаги распространения по каждому направлению и шаги затухания и вращения электрического поля сигнала. Также на шагах распространения по пространственным направлениям используется противопотоковая аппроксимация пространственных производных (метод Годунова с коррекцией потоков). Схема является консервативной, монотонной, имеет второй порядок точности по времени и третий по пространственным переменным.
Схема численного интегрирования уравнений Максвелла в литосфере
Пусть r =( x , y , z ) — декартовы координаты; t — время; E ( r , t ), D ( r , t ), H ( r , t ) и B ( r , t ) — напряженность и индукция электрического и магнитного полей; j ( r , t ) — плотность тока в момент времени t в точке с радиус-вектором г . Рассмотрим уравнение Фарадея и уравнение Максвелла в системе СИ
^ B = - rot E ( r , t ), ^ D = rot H ( r , t ) - j ( r , t ). (1)
о t d t
В литосфере эти уравнения замыкаются законом Ома в форме
j ( r , t ) = о ( г ) Е ( г , t ),
где o( r ) — скалярная проводимость среды, и материальными уравнениями
D ( r , t ) = S o S ( r ) E ( r , t ), B ( r , t ) = P 0 P ( r ) H ( r , t ), (3)
где e(r) и p(r) — безразмерные относительные ди- электрическая и магнитная проницаемости среды в низкочастотном пределе; е0 и р0 — электрическая и магнитная проницаемости вакуума.
В представленных расчетах мы полагали, что в литосфере p( r )=1, однако мы сохраним р в дальнейших формулах, поскольку используемая методика позволяет рассматривать общий случай с переменной p( r ).
Пусть c 0 = 1/ ^ £ 0 p 0 — скорость света в вакууме;
c(r) = c0 / V£(r)p(r) — скорость света в среде в точке с радиус-вектором r. Систему (1)-(3) можно записать в виде ав dE c2 +
= -rotE, = -rot д t д t £
B
Ц
- —E.
£ 0 £
Для численного интегрирования системы (4) применим метод расщепления по физическим процессам. Общий шаг интегрирования разделяется на два подшага. Один из них — это подшаг распространения, на котором интегрируется система уравнений dB д t
dE с22
=-rotE, = —rot II dt £(pj
второй — подшаг затухания сигнала, на котором аналитически интегрируется система уравнений dE _
--=--E dt£ по формулам E(r, t + t) = E(r, t) exp I - I. Пра-
( £0£(r)J вильное чередование подшагов расщепления обеспечивает второй порядок точности по времени.
Рассмотрим численное интегрирование системы
(5). Введем перенормированные поля
B(r, t) = c0B(r, t), E(r, t) = c0 / c(r)E(r, t), а также вектор M(r) = Vc в представленных расчетах, равный M(r) = Vc + cVц(r)/ ц(г) в общем случае. Умножая уравнение Фарадея на c0, получим уравнение
-
8 В
— = - rot(c E ). (6)
-
8 t
Подстановка введенных обозначений в уравнение Максвелла приводит это уравнение к виду dE = rot ( c — )-[m X — ].
Введем матрицы Rx, Ry, Rz и нулевую мат- рицу O:
|
" 0 0 |
0 ) |
f 0 0 1 ) |
||
|
R + = |
0 0 |
- 1 |
, R У = |
0 0 0 |
|
+ 0 1 |
0 у |
+- 1 0 0 , |
||
|
' 0 - 1 |
0' |
f 0 0 0 ) |
||
|
R z = |
1 0 |
0 |
, O = |
0 0 0 . |
|
+ 0 0 |
0 |
+ 0 0 0 v |
Оператор rot в декартовых координатах можно представить в виде
-
da da
rota = Rx — + Ry — + Rz —.
8x dy
Введем шестимерные вектор-столбец переменных u = ( B x , B y , B z , E x , E y , E z ) T и вектор F = ( F , F 2, F 3, F 4, F 5, F 6 ) T , компоненты которого заданы формулами
F j = F 2 = F 3 = 0, ( F 4 , F 5 , F 6 ) T = -[ m x — ] . (10)
Используя обозначения (8), а также формулы (9) и (10), систему уравнений (6), (7) можно представить в виде векторного уравнения du d(cu) d(cu) d( cu)
+ A x + A y + A z = F , (11)
8t 8x dy 8z где Ax, Ay, Az — симметричные матрицы размером 6X6, заданные формулами
^
A x
O R x
(- R x O )
^
О R y
(- R y O J
^
A z
'•> ^,
O R z
(- R z O J
Векторное уравнение (11) задает шестимерную линейную гиперболическую систему уравнений первого порядка, записанную в консервативной форме. Ее правая часть F линейно зависит от компонент вектора u. Для численного интегрирования таких систем разработано достаточно много разностных схем, в том числе схемы повышенного порядка точности, которые применяются для уравнений газовой динамики. Наиболее полное описание этих схем содержится в монографиях [Куликовский и др., 2012; Бисикало и др., 2013].
Используя метод расщепления по пространственным направлениям, можно построить явные монотонные схемы численного интегрирования системы (11), которые сводятся к последовательному интегрированию одномерных по пространству гиперболических систем уравнений. На каждом временном шаге нужно последовательно проинтегрировать три системы уравнений, например в таком порядке:
d u ' d ( c u )
+ A x— 8 t 8 x
= F , ^ A , ^ = f , x 8 t 8 y y
d um
---+
8 t
A z
d ( c "') = F d z z
При этом правые части систем (13) выбираются так, чтобы они не изменялись на своем шаге расщепления и удовлетворяли равенству
F = F x + F y + F z .
Эти условия будут выполнены, если задать их по формулам
T
0,0,0,0, - MB. , mb \ , zx yx
T
0,0,0, - M z B y ,0, - M x B y ) ,
T
0,0,0, - M y B z , M x B z ,0 ) .
На каждом шаге расщепления для двух компонент магнитного поля и двух компонент электрического поля, ортогональных направлению шага, рассчитывается распространение сигнала конечно-разностным способом. В качестве начальных условий для каждой системы уравнений в (13) берутся значения, рассчитанные на предыдущем шаге расщепления. Сохранить второй порядок аппроксимации по времени в схеме расщепления можно, если циклически изменять порядок выполнения шагов расщепления, например, выполняя сначала шаги по пространственным направлениям в следующем порядке: xyz , yxz , zxy , xzy , yzx , zyx . Обоснование этого утверждения содержится, например, в монографиях [Бисикало и др., 2013; Ковеня и др., 1981] .
Пусть заданы равномерная сетка по времени tn=10+nт, где т — шаг по времени, и равномерная пространственная сетка в декартовых координатах, целые и полуцелые узлы которых заданы соотношениями x = x0 + hxx, y, = y 0 + jhy , zk = z 0 + khz, r, j, k =( x,, y,, zk), xi+1/2 = x0 + (i + 1Z 2) hx, y,+1/2 = y0 + (j + 1Z 2) hy , zk+1/2 = z 0 +( k + 1/2) hz , в которых hx, hy, hz — шаги сетки по осям X, Y, Z соответственно. Будем использовать для значений функции f в узлах сетки обозначения f (rj,k, tn) = f ”,,k.
Рассмотрим послойный переход для первой системы из (13), которая интегрируется на шаге расщепления по направлению х . В этом случае представим заданную первой формулой в (10) матрицу
Ax в виде Ax = QЛQ , где Л — диагональная матрица, на диагонали которой стоят собственные числа матрицы Ax; Q — матрица, столбцы которой есть правые собственные векторы матрицы Ax, опре
- 1
деленные с точностью до множителя; Q — матрица, обратная Q. Строки матрицы Q 1 — левые соб ственные векторы матрицы Ax. Эти матрицы постоянны. Их можно взять в виде
( 0 0 0 0 0 0 )
0 1 0 0 0 0
Л =
0 - 1 0 0
0 0 0 0
0 0 0 0 1 0
( 0 0
0 0
0 - 1 )
(1 0 0 0 00
0 1/2 0 00 1/2
0 0 1/2 0 1/20
0 0 0 1 00
0 0 -1/2 0 1/20
( 0 - 1/2
0 0
0 1/2 )
0 0 0
0 0 ^
0 1 0 0
0 - 1
-•---1
Q
0 0
0 0
0 0
. 0 1
1 0 - 1 0
0 1
0 0
1 0 1 0
0 0 0
С помощью введенных матриц Л , Q и Q 1 первую систему уравнений в (13) представим в виде векторного уравнения
1 д ( c u ) д x
du '"
+ Q Л Q д t
= F x .
Введем вектор-столбец характеристических переменных (линеаризованные инварианты Римана)
1 по направлению x , заданный формулой w = Q u .
Компоненты вектора w заданы формулами
—-—~—•■ -"^—^ '•—~—^ .--——^ --—^—■ w1 Bx , w2 = By - Ez , w3 = Bz - Ey , w4 = Ex, w5 = Bz + Ey , w6 = By + Ez .
- 1
Умножим уравнение (15) слева на матрицу Q .
В результате получим векторное уравнение
^ w + Л ^w = Q - 1 F , (16)
д t дx x которое является системой из шести скалярных уравнений для характеристических переменных.
В итоге численное интегрирование первой системы из (13) сводится к численному интегрированию четырех независимых одномерных по пространству уравнений переноса (16) для характеристических переменных w2 = By - Ez , w3 = Bz - Ey , w5 = Bz + Ey, w6 = By + Ez .
Это интегрирование выполняется с помощью явной противопотоковой схемы, имеющей третий порядок точности по пространству и второй — по времени. Эта схема подробно описана в работах [Мингалев и др., 2018а, б] .
Рассмотрим кратко основные идеи, лежащие в основе численной схемы.
Схема численного интегрирования уравнений Максвелла в ионосфере
Для ионосферы мы полагали, что безразмерная относительная магнитная проницаемость среды p ( r ) = 1 и выполняется формула B ( r , t ) = ц 0 H ( r , t ). Также мы полагали, что имеет место поляризация плазмы D ( r , t ) = s 0 E ( r , t ) + P ( r , t ), где P ( r , t ) — вектор поляризации, причем плотность тока поляризации д Р / д t = j совпадает с полной плотностью тока в плазме. Уравнение Фарадея и уравнение Максвелла в этом случае принимают вид
— = - rot E , — = c 0'rot B - — j . (17)
д t , д t 0 s 0
Система (17) замыкается уравнением для плотности тока электронов, вызванного полем сигнала:
fj =- V ej -^e [ jX b^ro^E, д t где ve — частота столкновений электронов, Qe = e |Ввн| / me — гирочастота электронов, Ввн — внешнее геомагнитное поле, b = B / B , вн вн
® 2 = e 2 n e / ( m e e 0 ) — квадрат плазменной частоты электронов.
В ионосфере, как и в литосфере, применяется метод расщепления по физическим процессам. Правильное чередование подшагов расщепления обеспечивает второй порядок точности по времени. На подшаге распространения учитывается только распространение сигнала и численно интегрируется система уравнений
^ B = - rot E , — = c 22 rot B (18)
д t д t 0
с помощью той же схемы, что и система (5) в литосфере.
На подшаге затухания и вращения сигнала магнитное поле сигнала не изменяется и учитываются только затухание электрического поля и его вращение за счет внешнего геомагнитного поля. При этом в каждой точке расчетной сетки аналитически интегрируется система уравнений
| E = - — j, = - V e j a [ j X b ] + ro 2 S o E .
д t ^ o д t
Эта система разделяется на две независимых системы для продольных электрического поля E = ( b , j ) и тока j\ = ( b , j )
d E 2 . d j _ 2
, = j\\ , — = - V e j\\ + ® e E 0 E || (20)
dt Eo dt и поперечных электрического поля E± — E - bЕц и тока L = j — bj|| :
dE± = _2 • dt Eo
± ,
^ j T = - V e j ± - Q e [ L X b ] + ® 2 e 0 E ± .
d t
Системы (20) и (21) являются автономными линейными системами однородных дифференциальных уравнений (ОДУ) с постоянными коэффициентами. Решения задачи Коши для этих систем выражаются достаточно громоздкими аналитическими формулами.
В той части ионосферы, где выполнено неравенство ve t> 20, системы (20) и (21) заменяются на более простые уравнения и систему, дЕ\ д t
—
2 ® 2F
—Е\\
v e
5 E ± _ 2 ® 2 v e д t v 2 + Q 2
E ± — ^^ [ E ± x b ] v
которые также имеют аналитические решения задачи Коши.
Изложенная схема позволяет учесть частотную дисперсию и моделировать распространение широкополосных сигналов произвольной формы.
Условия на границе области моделирования
В модельных экспериментах на всех внешних границах было задано граничное условие свободного ухода волны, что достигалось обнулением входящего из-за пределов области моделирования потока. Такой подход позволяет получить низкие коэффициенты отражения плоской волны от границ области моделирования: для волн, падающих под углом от 80° до 90°, отношение амплитуды отраженной волны к амплитуде падающей плоской волны не превышает 0.01. При угле падения 60° это отношение составляет уже ~0.05, при угле 45° — ~0.16, при угле 27° — ~0.33, а при угле 18.4° — ~0.43 [Мингалев и др., 2018а]. Для сравнения: метод FDTD (finite-differences time-domain method) [Yee, 1966] при применении простых граничных условий, таких как условия Мура [Mur, 1981] и Лиао [Liao, 1984], дает отражения порядка 0.1…1 %, но только при падении волны на границу под прямым углом. При падении под острым углом коэффициент отражения растет до 100 % при падении по касательной. Однако при использовании непрерывно действующего источника даже столь малых отражений, которые порождает применяемая схема, достаточно для накопления ошибок в области моделирования, поэтому возникает необходимость использования методов подавления, подобных PML (perfectly matched layer), использующихся в FDTD-моделях [Berenger, 1994]. Именно такой тип источника применялся в представленных авторами экспериментах, что привело к необходимости адаптации и применения метода PML. Разделение схемы по пространственным переменным и физическим процессам позволяет применять профиль электрических и магнитных потерь, предложенный Беренгером, непосредственно к потокам противопотоковой схемы на границе области моделирования. Геометрический профиль потерь внутри отдельного слоя имеет вид
Р ( r) = — ^° c 0 ln( g )n ln( R o ) g* r / " ^ ) , (22) 2 A x ( g N — 1)
где g — коэффициент геометрической прогрессии; Δ x — шаг по пространству; c 0 — скорость света; N — номер PML-слоя, считая от интерфейса счетного региона и границы; r — расстояние от границы; R 0 — коэффициент отражения от первого слоя. В представленных численных экспериментах авторы используют профиль потерь, рассчитанный по формуле (22) со следующими параметрами: R 0 =0.01 (1 %), коэффициент прогрессии g =2.15, количество слоев 14. Несмотря на то, что коэффициент отражения от первого слоя не меньше, чем характерный для данной схемы при обнулении исходящих потоков на углах падения 80°–90°, а на практике даже больше вследствие отражений от последующих слоев, основным преимуществом метода PML является его крайне слабая зависимость от угла прихода электромагнитной волны. Данную особенность демонстрирует и адаптированный для противопотоковой схемы вариант.
РЕЗУЛЬТАТЫ И ОБСУЖДЕНИЕ
Авторами были проведены численные эксперименты семи типов: четыре эксперимента — с модельной областью, однородной по горизонтали и вертикальными профилями электронной плотности, как на рис. 1; три серии экспериментов — с переходом в направлении распространения волны, представленные на рис. 2. Модельный сигнал, исследуемый для всех представленных в работе численных экспериментов, задавался синусоидальной функцией с частотой 1500 Гц. Авторы не обсуждают переходные процессы, связанные с началом счета и прохождением фронта волны через области с не нулевым горизонтальным градиентом концентрации электронов. Все обсуждаемые далее результаты получены для устоявшегося режима (30-й –35-й периоды колебаний от начала моделирования). В связи с этим анализировать непосредственно временные формы компонент электромагнитного сигнала нецелесообразно, так как они идентичны исходному сигналу. Далее авторы сравнивают и обсуждают отдельные характеристики электромагнитного сигнала, такие как амплитуда и фаза, а также их производные.
Рис. 3. Разность амплитуд в процентах для электрического ( а ) и магнитного полей ( б ) для случаев распространения сигнала, когда вся область имела профиль электронной концентрации, характерный для спокойных условий 16.03.2013 в 06:00 UT и когда горизонтальный профиль концентрации электронов имел переход из спокойной области в область с повышенной концентрацией электронов
На рис. 3, а изображена разность амплитуд в процентах основной компоненты электрического поля волны для случаев распространения сигнала, когда вся область имела профиль электронной концентрации, характерный для спокойных условий 16.03.2013 в 06:00 UT и когда горизонтальный профиль концентрации электронов имел переход из спокойной области в область с повышенной концентрацией электронов, которая имела профиль, характерный для условий магнитной бури 17.03.2013 в 06:00 UT. На рис. 3, б изображена разность амплитуд основной компоненты магнитного поля волны в тех же экспериментах. На рис. 4, а , б изображена разность амплитуд основной компоненты электрического и магнитного полей волны для случаев распространения сигнала, когда вся область имела профиль электронной концентрации, характерный для условий магнитной бури 17.03.2013 в 06:00 UT и когда горизонтальный профиль концентрации электронов имел переход из области с повышенной концентрацией электронов в спокойную область, которая имела профиль, характерный для 16.03.2013, 06:00 UT. Видно, что распространение через области как положительного, так и отрицательного горизонтального градиентов электронной плотности оказывает влияние на амплитудные характеристики электромагнитных волн тем сильнее, чем выше градиент в области перехода. Однако в случае слабых возмущений электронной плотности изменения амплитуд основных компонент электромагнитного поля не достаточны для регистрации методами наземных наблюдений.
Анализ прохождения электромагнитной волной полосы повышенной концентрации электронов, расположенной перпендикулярно направлению распространения и параллельно плоскости волны, не выявил значительного влияния на амплитудные характеристики сигнала, независимо от соотношения ширины полосы и длины волны. Очевидно, что отдельная регистрация компонент электромагнитного поля не может быть использована для анализа слабых возмущений концентрации электронов в ионосфере, так как амплитудный отклик на них не превышает значений фонового шума.
Ситуация несколько улучшается в случае совместного анализа электрических и магнитных компонент электромагнитного поля. На рис. 5 показаны графики волнового сопротивления среды, рассчитываемого как отношение напряженностей магнитного и электрического полей для случаев горизонтально однородной среды с вертикальными профилями всех типов, представленных на рис. 1, и для случаев прохождения полосы повышенной концентрации электронов шириной λ и 0.25λ. Видно, что такой подход позволяет достаточно четко отделять случаи с наличием областей возмущений электронной плотности в ионосфере на радиотрассе от случаев спокойных условий, но слабо реагирует как на величину возмущений электронной плотности, так и на горизонтальную структуру данных возмущений. Для реализации данного подхода при анализе данных регистрирующая аппаратура в случае плоской волны, распространяющейся в плоскости волновода Земля—ионосфера, должна позволять фиксировать
Рис. 4. Разность амплитуд в процентах для случаев распространения сигнала, когда вся область имела профиль электронной концентрации условий магнитной бури 17.03.2013 в 06:00 UT и когда горизонтальный профиль концентрации электронов имел переход из области с повышенной концентрацией электронов в спокойную область: а — напряженность электрического поля; б — напряженность магнитного поля
Рис. 5. Волновое сопротивление среды
вертикальную компоненту электрического поля и две горизонтальные компоненты магнитного поля. В случае волн другого типа или произвольного направления распространения необходима фиксация всех шести компонент электромагнитного поля.
На рис. 6 показаны разности фаз напряженностей магнитного и электрического полей для случаев горизонтально однородной среды с вертикальными профилями всех типов, представленными на рис. 1, и для случаев прохождения полосы повышенной концентрации электронов шириной λ и 0.25λ. Видно, что фазовые характеристики наиболее сильно дифференцированы в зависимости от величины возмущения электронной плотности как в случае изменения вертикального профиля концентрации электронов, так и в случае горизонтального профиля в сравнении с ранее рассмотренными. Очевидно, что остается неопределенность, какого типа возмущение выз-
Рис. 6. Разность фаз между напряженностями электрического и магнитного полей
вало изменение фазы в каждом конкретном случае. Также остается вопрос зависимости сдвига фазы от формы профиля электронной концентрации, а также ее величины для волн различных диапазонов. Эти и многие другие вопросы требуют дополнительной проработки, выходящей за рамки представленной работы.
Для проведения анализа фазовых характеристик необходимо при регистрации компонент электромагнитного поля осуществлять точную привязку по времени регистрируемых электрических и магнитных компонент. В случае создания сети приемников, что может помочь с идентификацией горизонтальных ионосферных неоднородностей, необходимо обеспечить точную привязку по времени регистрируемых компонент поля к мировому времени, как, например, это делается в работе [Пильгаев и др., 2008] .
ЗАКЛЮЧЕНИЕ
Методами численного эксперимента авторами выполнена задача определения характеристик электромагнитных волн ИНЧ-диапазона, регистрируемых на уровне приземного слоя, несущих максимальное количество информации о состоянии волновода Земля—ионосфера.
Выявлено влияние горизонтальных градиентов концентрации электронов на амплитуду электромагнитных сигналов ИНЧ-диапазона при переходе между областями с различными характеристиками волновода Земля—ионосфера. Показано, что в случае незначительных изменений электронной плотности изменения амплитуды электромагнитного сигнала незначительны и не превышают 1 %.
Выявлено существенное влияние вида вертикального профиля электронной концентрации на волновое сопротивление среды и разность фаз между магнитной и электрической компонентами. Показано, что в случаях увеличения концентрации электронов в D-области волновое сопротивление среды увеличивается, а компонента электромагнитного поля E начинает запаздывать относительно H . По мнению авторов, отставание фазы электрического поля от магнитного обусловлено инерцией электронов и их столкновениями с нейтральными частицами ионосферы. Вызванная электрическим полем сигнала плотность электрического тока влияет на сдвиг фаз между моментами достижения максимальных значений электрическим и магнитным полями сигнала в фиксированной точке ионосферы.
На основании проведенных расчетов можно сделать вывод, что минимально пригодные для мониторинга состояния ионосферы наземные системы регистрации ИНЧ-волн должны содержать в своем составе три регистратора: один для вертикальной электрической компоненты и два для горизонтальных магнитных компонент электромагнитного поля. При этом важно обеспечить качественную синхронизацию между регистраторами и мировым временем.
Исследование выполнено при поддержке гранта Российского научного фонда (проект № 18-77-10018).
Список литературы Определение характеристик ИНЧ-волн, наиболее сильно реагирующих на незначительные изменения электронной плотности ионосферы в области высоких широт
- 1. Акимов В.Ф., Калинин Ю.К. Введение в проектирование ионосферных загоризонтных радиолокаторов. М.: Техносфера, 2017. 492 с.
- 2. Бисикало Д.В., Жилкин А.Г., Боярчук А.А. Газодинамика тесных двойных звезд. М.: Физматлит, 2013. 632 с.
- 3. Воскресенский Д.И., Степаненко В.И., Филиппов В.С. Устройства СВЧ и антенны. Проектирование фазированных антенных решеток. 4-е изд. М.: Радиотехника, 2012. 744 с.
- 4. Ковеня В.М., Яненко Н.Н. Метод расщепления в задачах газовой динамики. Новосибирск: Наука, 1981. 304 с.
- 5. Куликовский А.Г., Погорелов Н.В., Семенов А.Ю. Математические вопросы численного решения гиперболических систем уравнений. 2-е изд. М.: Физматлит, 2012. 656 с.
- 6. Мингалев И.В., Мингалев О.В., Ахметов О.И., Суворова З.В. Явная схема расщепления для уравнений Максвелла // Математическое моделирование. М.: РАН, 2018а. Т. 30, № 12. С. 17-38. DOI: 10.31857/S023408790001934-1.
- 7. Мингалев О.В., Мингалев И.В., Мельник М.Н., Ахметов О.И. и др. Новый метод численного интегрирования системы Власова-Максвелла // Математическое моделирование. М.: РАН. 2018б. Т. 30, № 10. С. 21-43. DOI: 10.31857/ S023408790001919-4.
- 8. Пильгаев С.В., Ахметов О.И., Филатов М.В., Федоренко Ю.В. Универсальное устройство синхронизации данных от GPS-приемника // Приборы и техника эксперимента. М.: РАН, 2008. № 3. С. 175-176.
- 9. Терещенко В.Д., Васильев Е.Б., Овчинников Н.А., Попов А.А. Средневолновый радиолокатор Полярного геофизического института для исследования нижней ионосферы // Техника и методика геофизического эксперимента, Апатиты: Изд. Кольского научного центра РАН, 2003. С. 37-46.
- 10. Berenger J.-P. A perfectly matched layer for the absorption of electromagnetic waves // J. Computational Phys. 1994. V. 114, N 2. P. 185-200. DOI: 10.1006/jcph.1994.1159.
- 11. Korja T., Engels M., Zhamaletdinov A.A., et al. Crustal conductivity in Fennoscandia - a compilation of a database on crustal conductance in the fennoscandian shield // Earth Planets Space. 2002. V. 54, N 5. P. 535-558. DOI: 10.1186/BF03353044.
- 12. Liao Z.P., Wong H.L., Yang B.P., Yuan Y.F. A transmitting boundary for transient wave analyses // Scienta Sinica (series A). 1984. V. 27, N 10. P. 1063-1076.
- 13. Mur G. Absorbing boundary conditions for the finite-difference approximation of the time domain electromagnetic field equations // IEEE Electromagnetic Compatibility. 1981. V. 23, N 4. P. 277-382. DOI: 10.1109/TEMC.1981.303970.
- 14. Yee Kane. Numerical solution of initial boundary value problems involving Maxwell’s equations in isotropic media // IEEE Transactions on Antennas and Propagation. 1966. V. 14. P. 302-307. DOI: 10.1109/TAP.1966.1138693.