DOI: 10.25702/^С.2307-5252.2018.9.5.157-164 УДК 519.6+537.877
О. И. Ахметов, И. В. Мингалев, О. В. Мингалев, З. В. Суворова
ВЛИЯНИЕ ПРОФИЛЯ ПРОВОДИМОСТИ ЛИТОСФЕРЫ НА РАСПРОСТРАНЕНИЕ ЭЛЕКТРОМАГНИТНЫХ СИГНАЛОВ СНЧ И КНЧ ДИАПАЗОНОВ В ВОЛНОВОДЕ ЗЕМЛЯ-ИОНОСФЕРА В ОБЛАСТИ ВЫСОКИХ ШИРОТ
Аннотация
Представлены результаты численных экспериментов прохождения электромагнитных сигналов СНЧ и КНЧ диапазонов в волноводе Земля -ионосфера в области высоких широт при различных профилях проводимости литосферы.
Ключевые слова:
численное моделирование, литосфера, СНЧ, КНЧ.
O. I. Akhmetov, I. V. Mingalev, O. V. Mingalev, Z. V. Suvorova
THE LITHOSPHERIC CONDUCTIVITY PROFILE INFLUENCE ON ELECTROMAGNETIC SLF AND ELF SIGNALS PROPAGATION IN THE HIGH LATITUDEEARTH - IONOSPHERE WAVEGUIDE
Abstract
The results of numerical experiments on propagation of the electromagnetic signals SLF and ELF ranges in the Earth-ionosphere waveguide in the region of high latitudes with different conductivity profiles of the lithosphere are presented.
Keywords:
numerical modeling, lithosphere, SLF, ELF.
Введение
Хорошо известно высокое значение ионосферных эффектов для всех видов деятельности, связанной с распространением радиоволн в волноводе Земля - ионосфера. Влиянием же второй границы волновода - литосферой, часто принято пренебрегать. Авторы представили в работе первые оценки эффективности данного подхода на Кольско-Карельском сегменте Балтийского щита сделанные методами численного моделирования распространения электромагнитных волн в волноводе Земля - ионосфера.
С целью выявления степени влияния литосферы были проведены эксперименты по прохождения электромагнитных сигналов в волноводе Земля -ионосфера при разных вариантах литосферных условий и фиксированных атмосферных. Вертикальный профиль концентрации электронов был одинаков для всех случаев, а атмосфера считалась однородной по горизонтали в области наблюдений, которая примерно соответствующей центральной части Кольского полуострова.
Описание численных экспериментов
Авторами были проведены расчёты распространения радиоволн в волноводе Земля-ионосфера для трех разных конфигураций литосферных условий. Атмосфера во всех случаях была идентичной. Вертикальный профиль концентрации электронов представлен на рис. 1а, он получен методами экстраполяции и интерполяции по данным модели 1Ы 2016 для даты 00:00иТ 28.08.2014. Во всех случаях модель была горизонтально однородной как в литосфере, так и в атмосфере, диэлектрическая проницаемость литосферы е принималась равной 9. В первом случае литосфера принималась однородной и по высоте с проводимостью а=110~5 См.
Во втором случае, представленном на рис. 1б, проводимость литосферы росла по экспоненциальному закону на глубинах более 10 км, а на глубинах менее 10 км принималась равной а=210~5 См. Данный профиль в целом соответствует результатам измерений, представленным в работе по исследованию проводимости грунта на Кольском полуострове [1].
Третий случай в целом идентичен второму за исключением того, что начальная проводимость приповерхностного слоя составляла а=210~3 См.
Источник во всех экспериментах представлял собой горизонтальный гармонический диполь длинной 100 км.
ю5 ю10 Ю"4 1°"3
Концентрация электронов м"3 Проводимость см
Рис. 3. Аппроксимированный высотный профиль концентрации электронов в 0 ИТ 28.08.2014 (а); аппроксимированный профиль проводимости литосферы (б).
Описание модели, граничных условий и функции источника
Используемая в представленной работе модель распространения электромагнитных сигналов в различных средах построена на основе схемы с противопотоковой аппроксимацией пространственных производных (метод Годунова с коррекцией потоков). Также используется расщепление по пространственным направлениям и по физическим процессам, причем затухание поля сигнала за счет проводимости и его вращение при наличии холловской проводимости среды учитываются на отдельных шагах расщепления по аналитическим формулам. Схема
является монотонной, имеет 2-й порядок точности по времени и 3-й по пространственным переменным, а также является консервативной.
Рассмотрим кратко основные идеи, лежащие в основе численной схемы,
более подробно в работах [2, 3].
% + ±Л1хи) + ±{1уи) + ±{1,и) = р (1)
где
и {Вх,Ву,В2,Ех,Еу,Е2), Р=(Р1,Р2,Р3,Р4,Р5,Р6У,
В = с0В, Е = фщЕ, ] =]
РХ = Р7=Р^ = 0,
Фи
Е0фр
"2
(Р4,Р5,Р6)Т = -[МхВ]-чЕ-],
М = Ус + 1
0=
■, с = С0/,17Л, Ч = а/(£0£)
£0 и - электрическая и магнитная постоянные, £ и ^ - безразмерные относительные диэлектрическая и магнитная проницаемости среды зависящие от пространственных координат. Симметричные матрицы Ах, Ау, Аг определяются формулами:
Ах =
О
сЛу
йхС
О
Ау =
О
сО
у
ИуС
О
О
сЯ,
0 0 0
0 0 1^
0
10
Ях = (0 0
1), Яу = ( 0 0 0), Ях = (1 0 0
0 1 0
100
000
Векторное уравнение (1) задает 6-мерную линейную гиперболическую систему уравнений 1-го порядка, записанную в консервативной форме. Ее правая часть ¥ линейно зависит от компонент вектора и. Для численного интегрирования таких систем разработано достаточно много разностных схем, в том числе схемы повышенного порядка точности, которые применяются для уравнений газовой динамики. Наиболее полное описание этих схем содержится в монографиях [4, 5]. Используя метод расщепления по пространственным направлениям и физическим процессам [5,6], можно построить явные монотонные схемы численного интегрирования системы (1), которые сводятся к последовательному интегрированию 1-мерных по пространству гиперболических систем уравнений.
Для изотропной среды на каждом временном шаге нужно последовательно проинтегрировать 3 системы уравнений, например, в таком порядке:
ди' д
д х
Ц- + д(Ауи") = Ру
д
ди'''
д у д
(2)
+ тЛА^и'") = Рг
д1 дг
При этом правые части систем (2) выбираются так, чтобы они не изменялись на своем шаге расщепления и удовлетворяли равенству Р = Рх + Ру + Рг .На каждом из шагов расщепления для двух компонент магнитного поля и двух компонент электрического поля, ортогональных направлению шага, рассчитывается распространение сигнала
и
(F4, F5, F6y = - (£) RXB - RyB - (£) RZB -c[MxB]-fjE-J
конечно-разностным способом, а также рассчитывается по аналитическим формулам затухание за счет проводимости третьей компоненты электрического поля. В качестве начальных условий для каждой системы уравнений в (2) берутся значения, рассчитанные в результате предыдущего шага расщепления. Сохранить второй порядок аппроксимации по времени в схеме расщепления можно, если циклически изменять порядок выполнения шагов расщепления. Например, выполняя сначала в следующем порядке шаги по пространственным направлениям: xyz, yxz, zxy, xzy, yzx, zyx. Обоснование этого утверждения содержится, например, в монографиях [5, 6].
В случае анизотропной среды с холловской проводимостью (такой средой является плазма в ионосфере и магнитосфере) тензор проводимости представляется в виде суммы его симметричной и антисимметричной частей. В этом случае к трем шагам расщепления схемы (2) добавляется четвертый шаг расщепления. Параметры среды и переменные, зависящие от проводимости и диэлектрической проницаемости, становятся тензорными. Компоненты вектора F задаются формулами:
г4, г5> F6)T = - (Ух) RXB - (д.y) RyB - (у) RZ
Т = £0-1£-1'2д£-1'2, c =
При этом на трех шагах расщепления (2) учитываются только симметричная часть тензора проводимости, а на четвертом шаге учитывается вращение электрического поля за счет антисимметричной части тензора проводимости, которое описывается аналитическими формулами. При этом магнитное поле не изменяется. На этом шаге в каждой точке расчетной сетки аналитически решается система уравнений:
§=[ПХЕ] (3)
здесь П = (fjyz, J]XZ, 7}ху) - угловая скорость, в которой fjyz, J]XZ, 7)xy являются компонентами антисимметричной части тензора Эта система задает вращение поля Е с вектором угловой скорости П . Соответствующее циклическое изменение последовательности выполнения шагов расщепления обеспечивает второй порядок аппроксимации по времени и в этом случае.
Во всех подставленных в работе численных экспериментах область моделирования представляла собой параллелепипед с основанием 1536x3072 км, высотой в атмосфере 300 км и глубиной в литосфере 100 км. Шаги сетки по горизонтали составляли 8 км, по вертикали в атмосфере 2 км и 1 км в литосфере. Шаг по времени составлял 4-10"6 сек.
В модельных экспериментах на всех внешних границах действовало граничное условие свободного ухода волны, что достигалось обнулением исходящего за приделы области потока. Применяемой авторами схема в таких условиях имеет достаточно низкие коэффициенты отражения плоской волны от границ области моделирования: для волн, падающих под углом от 80 до 90 градусов, отношение амплитуды отраженной волны к амплитуде падающей плоской волны не превышает 0.01. При угле падения 60 градусов это отношение уже составляет примерно 0.05, при угле падения 45 градусов - примерно 0.16, при угле падения 27 градусов - примерно 0.33, а при угле падения 18.4 градуса -примерно 0.43 [2]. Для примера метод FDTD (finite-differences time-domain method) [7] при применении простых граничных условий, таких как условия Мура (Mur) [8] и Лиао (Liao) [9] дает отражения порядка 0,1..1 %, но только при падении
волны на границу под прямым углом. При падении под острым углом коэффициент отражения растет вплоть до 100 % при падении по касательной. Однако при использовании непрерывно действующего источника даже столь малых отражений, которые порождает применяемая схема, достаточно для накопления ошибок в области моделирования, и возникает необходимость в использовании методов подавления подобных PML (perfectly matched layer), использующихся в FDTD моделях [10]. Именно такой тип источник применялся в представленных авторами экспериментах, что привело к необходимости адаптации и применения метода PML. Разделение схемы оп пространственным переменным и физическим процессам позволяет применять профиль электрических и магнитных потерь, предложенный Беренгером, непосредственно к потокам противопотоковой схемы на границе области моделирования. Геометрический профиль потерь внутри отдельного слоя имеет вид:
p(r) = -f^^Wo)g(r/Ax) (4)
где д - коэффициент геометрической прогрессии; Ах - шаг по пространству; с0 -скорость света; N - номер PML-слоя, считая от интерфейса счетного региона и границы; r - расстояние от границы; Ro - коэффициент отражения от первого слоя. В представленных численных экспериментах авторы используют профиль потерь, рассчитанный по формуле 4, со следующими параметрами: Ro=0.01 (1 %), коэффициент прогрессии д = 2.15, количество слоев 14. Несмотря на то, что коэффициент отражения от первого слоя не лучше чем характерный для данной схемы при обнулении исходящих потоков на углах падения 80-90 градусов, а на практике даже хуже вследствие отражений от последующих слоев, основным преимуществом метода PML является его крайне слабая зависимость от угла прихода электромагнитной волны. Данную особенность демонстрирует и адаптированный для противопотоковой схемы вариант.
В качестве источника сигнала во всех представленных экспериментах используется длинная линия (100 км), по которой протекает изменяющийся по закону синуса ток с максимальным значением 100 А. Линия направлена вдоль оси x, а по оси y смещена таким образом, чтобы расстояние до трех боковых границ счетной области от геометрического центра антенны было одинаково. Высота линии над поверхностью земли 15 м. Распределение поля в ближайших к линии 2-х узлах в каждую сторону рассчитывалось аналитически для условий £ = ^ = 1,а = 0.
Результаты и обсуждение
Результаты моделирования представляют собой массивы электромагнитных полей для всех моментов времени с модельным шагом и всех узлов сетки. Представлять результаты в таком виде в формате статьи не представляется возможным. Поэтому далее в обсуждении и на рисунках будут обсуждаться только огибающие сигналов полученные в результате усреднения максимумов амплитуды за два периода колебаний. Усреднение происходи спустя некий промежуток времени, различный для разных частот, когда начальные эффекты перестанут оказывать заметное воздействие на результирующие поля. Таким образом, от множества значений полей по времени мы перейдем к одному. По пространству результаты рассматриваются только в приземном слое в направлении перпендикулярном направлению антенны. Пример максимальных амплитуд электрического поля на четырех различных частотах в зависимости от расстояния до источника представлен
на рис. 2. Из рисунка видно, что с увеличением частоты сигнала скорость затухания с расстоянием увеличивается, структура поля вблизи источника отличается от структуры на значительном удалении. Такое поведение демонстрирует влияние волновода Земля - ионосфера на распространение электромагнитных сигналов.
Авторы проанализировали разности максимальных амплитуд полей для трех различных литосферных условий. Для первых двух случаев со сплошной литосферой проводимостью а=110~5 См и с профилем проводимости как на рис. 1 разница не превысила 1/1000 доли процента. Данный результат ясно демонстрирует, что в основном при проектировании антенн СНЧ, КНЧ диапазонов на Кольско-Карельском сегменте Балтийского щита влиянием структуры проводимости литосферы на прохождение сигнала в волноводе Земля-ионосфера допустимо пренебрегать и использовать для расчетов а=110~5 См. В случае значительно более высокой проводимости как в третьем модельном варианте литосферы а=210~3 См поведение сигнала изменяется как в зависимости от расстояния до источника, так и в зависимости от частоты. На рис. 3 в разных панелях представлены отклонения максимальной амплитуды электрического поля в зависимости от расстояния на разных частотах, для случаев два (профиль на рис. 1) и три (начальная а=210~3 См). Из рис. 3 для частот 5 и 10 Гц видно, что вблизи источника повышенная проводимость литосферы приводит к ослаблению поля на 0.5-1%, в то же время на значительном расстоянии более 400 км в случае 5 Гц и 300 км в случае 10 Гц, приводит к улучшению отражения от поверхности Земли и увеличению амплитуды. Для более высоких частот заметна обратная зависимость для компоненты Ег. Более значительные вариации компоненты Ех связаны с тем, что и сама эта компонента больше см. рис. 2. В целом с ростом частоты влияние профиля проводимости литосферы снижается.
Расстояние км Расстояние км
Рис. 2. Максимальные амплитуды компонент электрического поля на четырех различных частотах в зависимости от расстояния до источника.
Рис. 3. Отклонения максимальной амплитуды компонент электрического поля в зависимости от расстояния на разных частотах.
Заключение
В работе показано, что в основном при проектировании антенн СНЧ, КНЧ диапазонов на Кольско-Карельском сегменте Балтийского щита влиянием структуры проводимости литосферы на прохождение сигнала в волноводе Земля-ионосфера допустимо пренебрегать и использовать для расчетов а=1-10"5 См, кроме особенных случаев, к примеру, на побережье, где значительно более высокая проводимость литосферы. В этих случаях поведение сигнала изменяется как в зависимости от расстояния до источника, так и в зависимости от частоты.
Литература
1. Korja T., Engels M., Zhamaletdinov A.A., Kovtun A.A., Palshin N.A., Smimov M.Yu., TokarevA.D., AsmingV.E., VanyanL.L., VardaniantsI.L., and the BEAR Working Group. Crustal conductivity in Fennoscandia—a compilation of a database on crustal conductance in the Fennoscandian shield // Earth Planets Space. 54.2002.535-558.
2. Мингалев И.В., Мингалев О.В., Ахметов О.И., Суворова З.В. Явная схема расщепления для уравнений максвелла // Математическое моделирование. 2018. том 30. № 12. с. 17-38.
3. Мингалев О.В., Мингалев И.В., Мельник М.Н., Ахметов О.И., Суворова З.В. Новый метод численного интегрирования системы Власова-Максвелла // Математическое моделирование. 2018. том 30. № 10. с. 21-43.
4. Куликовский А.Г., Погорелов Н.В., Семенов А.Ю. Математические вопросы численного решения гиперболических систем уравнений. 2-е изд., испр. и доп. - М.: ФИЗМАТЛИТ, 2012. 656 с.
5. Бисикало Д.В., Жилкин А.Г., Боярчук А.А. Газодинамика тесных двойных звезд. Москва: ФИЗМАТЛИТ, 2013. 632 с.
6. Ковеня В.М., Яненко Н. Н. Метод расщепления в задачах газовой динамики. Новосибирск: Наука, 1981.
7. Yee Kane. Numerical solution of initial boundary value problems involving Maxwell's equations in isotropic media // IEEE Transactions on Antennas and Propagation. 1966. Vol. 14. PP. 302-307.
8. Mur G. "Absorbing boundary conditions for the finite-difference approximation of the time domain electromagnetic field equations", IEEE Trans. Electromagnetic Compatibility, vol. 23, no. 4, pp. 277-382, Nov. 1981.
9. Liao Z.P., Wong H.L., Yang B.P., Yuan Y.F. A transmitting boundary for transient wave analyses. In: Scienta Sinica (series A). 17.1984. S. 1063 - 1076
10. Berenger J-P. A perfectly matched layer for the absorption of electromagnetic waves // Journal of Computational Physics. Vol. 114. Is. 2. P. 185-200.
Сведения об авторах Ахметов Олег Иршатович
к.ф.-м.н., н. с., Полярный геофизический институт, Апатиты E-mail: akhmetov@pgia.ru
Мингалев Игорь Викторович
д.ф.-м.н., в. н. с., Полярный геофизический институт, Апатиты; E-mail: mingalev_i@pgia.ru
Мингалев Олег Викторович
к. ф.-м. н., зав. сектора, Полярный геофизический институт, Апатиты E-mail: mingalev_o@pgia.ru
Суворова Зоя Викторовна
програмист, Полярный геофизический институт, Апатиты E-mail: suvorova@pgia.ru