CC BY
2
Вестник КРАУНЦ. Физ.-мат. науки. 2023. Т. 44. №3. C. 104-120. ISSN 2079-6641
МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ d https://doi.org/10.26117/2079-6641-2023-44-3-104-120 Научная статья
Полный текст на русском языке УДК 537.874, 519.6
Компьютерная симуляция распространения плоской электромагнитной волны в волноводе образованном поверхностью Земли и ионосферой при условии неоднородной
проводимости границ
Д. А. Твёрдый*, Е. И. Малкин
Институт космофизических исследований и распространения радиоволн ДВО РАН, 684034, c. Паратунка, ул.Мирная, д. 7., Россия
Аннотация. В статье приводятся результаты компьютерного моделирования распространения плоской электромагнитной волны. Грозовые разряды являются естественными источниками импульсного электромагнитного излучения (атмосферика). Атмосферик распространяется в волноводе Земля-Ионосфера как плоская электромагнитная волна с широким спектром частот с максимумом интенсивности в диапазоне 4-10 кГц. После землетрясений с магнитудами порядка 7 и более, увеличивается насыщенность подземных вод вторичными минералами, что приводит к локальному увеличению проводимости земли в данных областях. Что определяет электрофизические свойства земли, как параметры нижней границы волновода Земля-Ионосфера. Которые влияют на характеристики распространяющихся в волноводе электромагнитных волн. Предполагается, что, изучая параметры атмосферика можно установить наличие неоднородности проводимости стенки волновода. На основе системы уравнений Максвелла с краевыми условиями в виде Perfectly matched layer, задаётся математическая модель процесса. Краевые условия модели определяют область распространения электромагнитного широкополосного сигнала, как волновод с неоднородной проводимостью границ. Система модельных уравнений решается численным методом Finite-Difference Time-domain. Для решения поставленной задачи и проведения компьютерного моделирования разработан программный комплекс в среде MATLAB. С целью проверки предположения, был проведён ряд компьютерных симуляций. В результате показано, что имеет место обратное рассеяние электромагнитной волны на трассе волновода, возникающее в следствии отражения волны при её взаимодействии с неоднородностью проводимости нижней границы волновода. Показано, что с помощью математического моделирования процесса распространения атмосферика и его взаимодействий с неоднородностью в волноводе можно установить наличие неоднородности и её связь с характеристиками излучения.
Ключевые слова: математическое моделирование, динамические процессы, уравнения Максвелла, атмосферик, плоская ЭМ волна, волновод, неоднородность проводимости, компьютерная симуляция, MATLAB, обратная волна, численное решение, FDTD, PML, ABC
Получение: 19.10.2023; Исправление: 27.10.2023; Принятие: 29.10.2023; Публикация онлайн: 02.11.2023
Для цитирования. Твёрдый Д. А., Малкин Е. И. Компьютерная симуляция распространения плоской электромагнитной волны в волноводе образованном поверхностью Земли и ионосферой при условии неоднородной проводимости границ // Вестник КРАУНЦ. Физ.-мат. науки. 2023. Т. 44. № 3. C. 104-120. EDN: XZYCQD. https://doi.org/10.26117/2079-6641-2023-44-3-104-120.
Финансирование. Работа выполнена в рамках Государственного задания по теме (2021—2023 гг.) "Физические процессы в системе ближнего космоса и геосфер при солнечных и литосферных воздействиях регистрационный номер АААА-А21-121011290003-0
Конкурирующие интересы. Конфликтов интересов в отношении авторства и публикации нет. Авторский вклад и ответственность. Авторы участвовали в написании статьи и полностью несут ответственность за предоставление окончательной версии статьи в печать.
* Корреспонденция: А E-mail: tverdyi@ikir.ru
Контент публикуется на условиях Creative Commons Attribution 4.0 International License © Твёрдый Д. А., Малкин Е. И., 2023
© ИКИР ДВО РАН, 2023 (оригинал-макет, дизайн, составление)
Vestnik KRAUNG. Fiz.-Mat. nauki. 2023. vol. 44. no. 3. P. 104-120. ISSN 2079-6641
MATHEMATICAL MODELING
" https://doi.org/10.26117/2079-6641-2023-44-3-104-120 Research Article Full text in Russian MSC 78A40, 93A30, 65M06
Computer Simulation of the Propagation of a Plane Electromagnetic Wave in a Waveguide Formed by the Earth's Surface and the Ionosphere Under the Condition of Inhomogeneous Boundary Conductivity
D. A. Tverdyi*, E.I. Malkin
Institute for Cosmophysical Research and Radio Wave Propagation FEB RAS, 684034, Paratunka, Mirnaya st., 7., Russia
Abstract. The article presents the results of computer simulation of the propagation of a plane electromagnetic wave. Lightning discharges are natural sources of pulsed electromagnetic radiation (atmosphere). Atmospheric propagates in the Earth-Ionosphere waveguide as a plane electromagnetic wave with a wide frequency spectrum with an intensity maximum in the range of 4-10 kHz. After earthquakes with magnitudes of the order of 7 or more, the saturation of groundwater with secondary minerals increases, which leads to a local increase in the conductivity of the earth in these areas. What determines the electrophysical properties of the earth, as the parameters of the lower boundary of the Earth-Ionosphere waveguide. Which affect the characteristics of electromagnetic waves propagating in the waveguide. It is assumed that, by studying the parameters of the atmosferic, it is possible to establish the presence of an inhomogeneity in the conductivity of the waveguide wall. Based on the system of Maxwell equations with boundary conditions in the form of a Perfectly matched layer, a mathematical model of the process is set. The boundary conditions of the model determine the region of propagation of an electromagnetic broadband signal as a waveguide with inhomogeneous boundary conductivity. The system of model equations is solved by the numerical Finite-Difference Time-domain method. In order to solve the problem and conduct computer modeling, a software package was developed in MATLAB environment. In order to verify the assumption, a number of computer simulations were carried out. As a result, it is shown that there is a backscattering of electromagnetic wave on the waveguide trace, arising as a consequence of reflection of the wave in its interaction with the inhomogeneity of the conductivity of the lower boundary of the waveguide. It is shown that with the help of mathematical modeling of the process of atmospheric propagation and its interactions with inhomogeneity in the waveguide it is possible to establish the presence of inhomogeneity and its relationship to the characteristics of radiation.
Key words: mathematical modeling, dynamic processes, Maxwell's equations, atmospheric, plane EM wave, waveguide, conduction inhomogeneity, computer simulation, MATLAB, inverse wave, numerical solution, FDTD, PML, ABC
Received: 19.10.2023; Revised: 27.10.2023; Accepted: 29.10.2023; First online: 02.11.2023
For citation. Tverdyi D. A., Malkin E.I. Computer simulation of the propagation of a plane electromagnetic wave in a waveguide formed by the Earth's surface and the ionosphere under the condition of inhomogeneous boundary conductivity. Vestnik KRAUNC. Fiz.-mat. nauki. 2023, 44: 3,104-120. EDN: XZYCQD. https://doi.org/10.26117/2079-6641-2023-44-3-104-120.
Funding. The work was carried out within the framework of the State task on the topic (2021—2023) "Physical processes in the system of near space and geospheres under solar and lithospheric impact registration number AAAA-A21-121011290003-0
Competing interests. There are no conflicts of interest regarding authorship and publication.
Contribution and Responsibility. All authors contributed to this article. Authors are solely responsible for providing the final version of the article in print. The final version of the manuscript was approved by all authors.
* Correspondence: A E-mail: tverdyi@ikir.ru ^jj
The content is published under the terms of the Creative Commons Attribution 4.0 International License © Tverdyi D. A., Malkin E.I., 2023
© Institute of Cosmophysical Research and Radio Wave Propagation, 2023 (original layout, design, compilation)
Введение
По данным WWLLN (The World Wide Lightning Location Network) каждую секунду на земле возникает до 50 молний. Молниевые разряды в грозовых облаках генерируют электромагнитное излучение естественного происхождения. Это позволяет рассматривать молниевые разряды как постоянно действующие источники импульсного электромагнитного излучения (атмосфериков) [1-3]. Атмосферик является широкополосным сигналом с максимумом интенсивности в диапазоне 6-10 кГц, что позволяет отнести их к классу сигналов очень низких частот (VLF) [4]. Однако, атмосферик распространяется не в свободном пространстве, а в сложной проводящей структуре волновода [5].
Как известно, волновод это структура направляющая электромагнитные или звуковые волны, с минимальными потерями энергии, ограничивая их передачу в одном направлении [6, 7]. В свою очередь, ионосферный волновод это пространства между слоями ионосферы, образующийся благодаря отражению радиоволн определённых частот от ионосферных слоёв, в которых плотность плазмы немонотонно меняется по высоте, образуя локальные максимумы и минимумы, за счёт чего происходит локализованное распространение радиоволн [8,9].
Рис. 1. Геометрия распространения лучей в ионосферном волноводе (www.soest.hawaii)
[Figure 1. Geometry of Ray Propagation in an Ionospheric Waveguide
(www.soest.hawaii) ]
В приземном ионосферном волноводе как на рисунке 1, который частично моделируется в данном исследовании, нижней границей служит поверхность Земли, а верхней - один из слоёв ионосферы [8].
При распространении атмосферика [10, 11] происходит неизбежное взаимодействие с неоднородностями в волноводе. Как правило, в большинстве научных исследований связанных с моделированием и изучением электромагнитной динамики [12-14] и многих других, рассматриваются именно неоднородности верхней стенки волновода т.е одного из слоёв ионосферы. В
данном исследовании нас интересует влияние неоднородностей именно нижней границы на распространение сигнала т.е. изменения проводимости на земной поверхности, которые существовать постоянно (прибрежная линия океанов), или быть временно возникающими (локальное изменение проводимости).
Последние, могут быть вызваны изменением химического состава подземных приповерхностных вод [15]. Известно что вода хороший проводник. Однако, если в составе вод возможно появление вторичных минералов в следствии геологической активности, или растворённых в водах свободных газов, например неэлектропроводного С02 [16]. То такого рода примеси в составе подземных приповерхностных вод могут быть причиной возникновения неоднородностей проводимости. Это может помочь связать некоторые гидрогеохимические предвестники землетрясений за недели-месяцы до их возникновения, с симулируемыми неоднородностями земли как нижней границы волновода и вызываемыми ими особенностями распространения атмосферика.
Математический аппарат моделирования
Уравнения Максвелла - это четыре фундаментальных математических выражения для описания электромагнитного поля, разделяемого на Е - вектор напряженности электрического поля в [вольт/метр] и Н - вектор напряженности магнитного поля в [ампер/метр]. Как правило, если необходимо провести моделирование распространения ЕМ волн то обращаются именно к этим известным уравнениям [17, 18], т.к. они связывают Е и Н поля корректно описывая их распространения при взаимном воздействии друг на друга. В дифференциальной форме записи, уравнения имеют вид:
Ух 7 = ™ + 7, (1а)
Ух ?? = -д? - М, (1Ь)
У- - = Ре, (1С)
V-B = pm. (1d)
Далее, для решения модельной задачи на основе уравнений Максвелла (1), мы будем использовать Finite-Difference Time-Domain (FDTD) численный метод на равномерных пространственных сетках [19]. Этот метод был одним из первых для численного решения (1) на пространственных сетках и оставался предметом постоянного развития. Метод на основан по "шахматному"принципу т.к. позволяет попеременно вычислять и использовать как E так и H значения полей, что делает решение более надежным, чем при использовании только одного из них.
Замечание 1. Тот факт, что уравнения Максвелла (1) на самом деле являются гиперболическими PDE, является основой всего метода FDTD. Т.е. для (1) есть решения которые ведут себя как волны, а возмущения распространяются с конечной скоростью.
За более чем 50 лет развития и применения, FDTD зарекомендовал себя с самой лучшей стороны, во многих задачах вычислительной электромагненики, в том числе и исследованиях близких к нашему. Например в исследовании Wang Z., Zhou C., и др. [20] применяется алгоритм Finite-Difference Time-Domain (FDTD) для моделирования глобального процесса распространения сейсмической EM волны в атмосферу и полость ионосферного волновода между нижней частью ионосферы и поверхностью земли. Предложенная модель распространения EM ELF связи литосфера-атмосфера-ионосфера очень перспективна для объяснения аномальных возмущений результатов наземных и космических наблюдений электромагнитных полей ELF, связанных с землетрясениями.
Замечание 1. Важной особенностью численного метода FDTD в том, что вычислительная сетка является бездивергентной по своим полям E и H в отсутствии свободных электрических и магнитных зарядов. Т.е. далее мы можем учитывать только уравнения ротора (1a-1b) при выводе схем FDTD, потому что уравнениям дивергенции (1c-1d) уже будут неявно учитываться в разработанных схемах FDTD [21, p.60].
Плотность электрического тока J из (1a) является суммой плотности тока проводимости J c = <ce E и плотности приложенного тока J i при J = J c + J ^ Аналогично в (1b), для плотности магнитного тока M = Mc + M^ где Mc = cmH.
Замечание 2. Следует отметить такие параметры физической среды, как:
• ce — электропроводность в [сименс/метр];
• cm — магнитная проводимость в [ом/метр].
так как с их помощью будут определяться границы проводящего объёма волновода Земля-ионосфера [7].
Для моделирования распространения сигнала, как изменения компонент полей E и H, уравнения (1) дополняются определяющими соотношениями D = e E а так же B = цН. Где e - диэлектрическая проницаемость [фарад/метр], ц — магнитная проницаемость материала [генри/метр]. Причём в свободном пространстве: e = e0 ~ 8.854 х Ю-12, ц = ц = 4п х Ю-7.
В итоге, теперь учитывая для (1a-1b) определяющие соотношения и правила разложения токов, мы можем представить уравнения в виде:
_^ d E _^ _^
V х H = e— + ce E + J i, (2a)
ot
V х lE = -ц^Д - cm H - Mi, (2b)
ot
в такой более упрощенной формулировке(2), присутствуют все четыре определяющих параметра е, ц, сте и стт, а значит можно задать любой линейный изотропный материал для моделирования распространения излучения в нём [22, р.3].
Метод РЭТЭ, так же называемый алгоритмом Йи [19], делит геометрию задачи на пространственную сетку, что позволяет разместить компоненты Е и Н полей в дискретных положениях пространства, и соответственно решать уравнения (2) в дискретные моменты времени. Поэтому далее, необходимо разложить каждое из двух уравнений (2) на три скалярных уравнения в декартовой системе координат для трёхмерного пространства.
Замечание 3. Соответственно поле Е будет разложено на Ех, Еу и Е2 компоненты, а Н на Нх, Ну и Н компоненты по направлениям в пространстве. Физические параметры среды е, ц - проницаемости, сте - электро и стт - магнитной проводимости материала, задаются непосредственно по каждому из направлений (Х,у,г) пространства, как СТреХ, ^рт^, ^рт%.
Из которых можно скомбинировать, в свою очередь две различные подзадачи: ТЕ и TM - моды, т.е. поляризации волн, по каждому из направлений для осей (х,у,г). Подробнее с выводом уравнений из (2) можно ознакомиться в работе [23].
О граничных условиях для задания области моделирования как волновода
Задание граничных условий позволят определить конкретную и конечную область моделирования. Последний пункт также важен, т.к. из-за особенности численной схемы FDTD [22, p.24] у такой сетки нет края, и его необходимо задать с помощью определенных условий. Например можно задать для компонент полей граничные значения равные 0. Но это создаёт свои особенности, а именно, такие граничные точки где значения E - поля, всегда равны нулю, начинают имитировать идеальный электрический проводник т.е. PEC (Perfect Electric Conductor), полностью отражая падающие волны обратно в пространство моделирования [24].
Однако, существуют условия усечения расчетной сетки призванные поглощать распространяющиеся наружу значения компонент полей с незначительным отражением, называемые ABC (Absorbed Boundary Condition). Для нашей задачи мы будем задавать ABC с помощью PML (Perfect matched layer) поглощающих сред конечной толщины, основанную на ложных определяющих параметрах для создания условия согласования волновых импедансов cre и crm, подробнее в [21, p.273]. PML располагается на краях вычислительного пространства, определённого вычислительной сеткой, подобно стенкам безэховой камеры [25].
Замечание 4. Мы будем использовать для задания неоднородности на нижней (Земля) стенке волновода Земля-ионосфера комбинацию этих двух способов усечения:
• PML будет имитировать распространение большей части излучения за границу волновода вдоль всего периметра моделируемой области;
• область неоднородности будет задаваться также PML областью со сложной геометрией и определёнными свойствами, позволяющими чуть большей части излучения отразиться от PEC стенок обратно на трассу волновода. Тот факт, что уравнения Максвелла (1) на самом деле являются гиперболическими уравнениями в частных производных, является основой всего метода FDTD. Т.е. для (1) есть решения которые ведут себя как волны, а возмущения распространяются с конечной скоростью [26].
ABC типа PML достаточно надёжны и эффективны в задачах вычислительной электромагнетики [27-29]. А также хорошо подходят для усечения решеток FDTD в свободном пространстве, в материалах с потерями или дисперсией, а также в металлических или диэлектрических волноводах. Однако, согласно [25] при задании постоянного однородного а по всей среде PML, мы получим заметные отражения подобные PEC, когда фронт электромагнитной волны в компьютерной симуляции подойдет к границе «PML - проводящий объём».
Поэтому, сопротивления ае и ат в среде PML следует задавать в виде распределения по некой функции, от малых значений около 0 на границе «PML -проводящий объём» до некоего заданного большего у противоположной стенки PML среды. Для вычисления такого профиля проводимости а(р) определим степенной закон распределения согласно [22,27]:
r , (- Y'" (npml + 1)£0Cln (R(0)) а(р) = O-max^J , amax =--2ASN-' (3)
где, R(0) - коэффициент отражения среды PML конечной ширины при нормальном падении; N - количество ячеек определяющих ширину PML; As - размер ячейки PML; npml - полиномиальный порядок оценки а(р), при npml = 1 - линейная, а при npml = 2 - параболическая; -Р - профиль потерь.
Значения магнитных проводимостей можно легко получить из отношений apmz = —apez и apmx = —apex, согласно [22, p.188]. Для am распределения по
профилям потерь аналогичны геометрически, на интервале [0,0.17].
Неоднородность на (рис. 2 a) определяется аналогично по формуле (3) при npml = 2, и N = 120, но уже при R(0) = 2e-1, а также при помощи боле сложных алгоритмов для указания области неоднородности и её возможной сложной геометрии.
Замечание 5. Тем самым, в данной математической модели, при R(0) = 2e-1 в такой зоне неоднородности будет иметь место пониженная сопротивляемость, что тождественно равно повышенной проводимости. В свою очередь, такие области повышенной проводимости имеют физические предпосылки [15] к возникновению в реальной среде.
а) Ь)
10 6 у 1(1 6
ПХ 0 О nz ПХ Ü О
Рис. 2. Профиль проводимости cre для неоднородности 1-м приближении. Пример распределения проводимостей по всей вычислительной сетке, включая распределения для PML по формуле (3) при R(0) = 1e-8, npml = 2, и N = 120: a) cre вдоль оси x; b) cre вдоль оси z [Figure 2. Conductivity profile cre for inhomogeneity of the 1st approximation. An example of the distribution of conductivities over the entire computational grid, including distributions for PML according to the formula (3) with R(0) = 1e-8, npml = 2, and N = 120: a) cre along the x axis; b) cre along the z axis ]
Единая численная схема описывающая распространение сигнала внутри волновода и его ГЫЬ стенках
Далее при выводе численной схемы мы будем отталкиваться от ТЕУ моды (с Ег, Ех и Ну компонентами полей) для моделирования распространения в 2-Б пространстве компоненты Ег напряженности электрического поля, которая является одной из характеристик электромагнитного изучения в волноводе.
Для удобства моделирования, классически область РМЬ считается является анизотропной средой или свободным пространством, полагая что во всёй области РМЬ параметры е = е0, Ц = Ц0. Также для упрощения отбросим поля независимых источников энергии, полагаем | 8оитсе = 0 и М5оитсе = 0.
Подробнее с выводом схемы для расчётов в области РМЬ из (2) можно ознакомиться в работе [23]. Но стоит отметить, что из-за моделирования ТЕУ -моды в 2-мерной плоскости, для вычисления распространения сигнала в области РМЬ, возникает необходимость в разбиении компоненты Ну на две Ну = Нуг + Нух искусственные, которые будут связанны с направлениями г и х. Компоненты поля Ех и Нуг вместе могут представлять волну, распространяющуюся в направлении г, а компоненты поля Ег и Нух представляют волну, распространяющуюся в направлении х. А так же то, что для среды РМЬ важны лишь значения проводимостей ае и ат, а именно, что РМЬ в отношении параметров е = е0, ц = ц0 - свободное пространство.
В случае такой постановки задачи, из (2) можно вывести общие уравнения для расчёта значений компонент Ег, Ех и Ну, и для внутреннего проводящего объёма и РМЬ границ волновода, которые примут вид. Тогда скалярные уравнения по
(х,у,г) осям для ТЕу - моды, примут вид:
9Ez + _ _ 9 (HyZ + HyX)
£ ~ , + opex Ez — ~ ,
dt 9x
9Ex + E __ 9 (Hyz + Hyx) £ 9t +0pezEc _ 9z (4)
ЭИуг + н _ — дЕ Ц + _ 9z,
ЭИух + И _ 9Е^ Ц 91 + артхИух_ ах"
Далее, данные уравнения можно представить как численные схемы РЭТЭ для расчёта компонент полей напряженности электрического Е и магнитного И полей:
En+'(i,k)_ Ceze(i,k) X E£(i,k) + Cezhy(i,k) x (Hn+1/2(i,k)_ H^2 (i,k _ 1)
+Hx+1/2(i,k)_ H^x+1/2(i,k _ 1)), (5)
2£z _ Atffpex(i,k) ^ (i лЛ 2At
Ceze(i, k) — -« л , л (^ЬгД'Цk) —
2£z + AtCTpeX (i,k)' zy ' (2£z + AtCTpex(i,k)) Ax'
En+1(i, k() _ Cexe(i,k) X En(i,k) +
Cexhy (i, k) X (Hn+1/2(i,k) _ Hn+1/2(i_ 1,k)
+Hnx+1/2(i,k)_ H^x+1/2(i _ 1,k)), (6)
2£x _ Atapez(i,k) плл_ 2At
Cexe(i, k) — -« л f. -i , (-exb^^ k) —
2£x + Atffpez(i,k^ y (2£x + Ato-pez(i,k)) Az
Hn+1/2(i,k) _ Chyzh(i,k) X Hn_1/2(i,k)+ Chyzex(i,k) X (En(i + 1,k) _ En(i,k)), (7)
_ AtOpmz(i, k) 2At
^-ЬХ^Ь^'Ц k) -« д, r. \ , Chг|zex(i, k)
+ AtOpmz(i, k)' yzex ' (2^y + AtOpmz(i, k)) Az'
H^x+1/2(i,k) _ Chyxh(i,k) x H^x_1/2(i,k)+ Chyxez(i,k) X (En(i,k + 1)_ En(i,k)), (g)
_ AtOpmx(i,k) а лл _ 2At
Chyxh(i, k) — -« л f. -i , Chyxez(i, k) —
+ А1арШх(г, к)' ухег ' (2цу + Агстртх(г, к)) Ах'
где Ах, Аг - приращения пространства решетки в направлениях координат х,г, соответственно, целые числа обозначающие номер ячейки; А1 - приращение времени, которое предполагается равномерным на интервале наблюдения; п — целое число; Индексы (г, к) - это размерность матриц, определяющих значения свойств физической среды и компонентов полей Е и И в пространстве 2-мерной задачи.
Параметры компьютерных симуляции
Для математического моделирования процесса взаимодействия электромагнитного излучения с неоднородностью в волноводе разработан программный комплекс в среде МЛТЬЛВ, реализующий описанный численный метод РЭТЭ и визуализирующая результат. В данной программе задаётся моделируемый волновод Земля-Ионосфера согласно схеме, как на рисунке 3 ниже, а значения полей Е и Н в каждой точке данной области в каждый момент времени вычисляется по формулам (5) - (8).
Рис. 3. Схематичное представление волновода Земля-ионосфера как вычислительной сетки [Figure 3. Schematic representation of the Earth-ionosphere waveguide as a
computational grid ]
где пх = 600 [ячеек] - размер моделируемого волновода вдоль Ох; пг = 2400 [ячеек] - размер моделируемого волновода вдоль Oz; мера длины [] - как число шагов дискретизации; n_pml_zn = 60 и n_pml_zp = 60 [ячеек] - толщина слоёв РМЬ имитирующего распространение широкополосного сигнала за область моделируемого волновода без отражений; п_рт1_хр = 60 [ячеек] - толщина слоя РМЬ определяющего верхнюю границу (Ионосфера) волновода Земля-Ионосфера; п_рт1_хп = 60 [ячеек] - толщина слоя РМЬ определяющего нижнюю границу (Земля) волновода Земля-Ионосфера; п_рш1_хп_Б = 240 [ячеек] - ширина области неоднородности нижней границы волновода вдоль Ог; , Б2("), 53("Ь) -
положения источника (9) импульсного электромагнитного излучения.
Для моделирования естественного источника импульсного электромагнитного излучения (атмосферика) будем использовать Гауссовский импульс, функция которого имеет вид:
( (" - "0\Л N0^23тах(Аг,Ах)
^«^«РЫ — ) , Т =-—-, (9)
где T = 2400 - общее время моделирования; t е [0, T] - текущее время моделирования; t0 - фаза; pos - номер позиции источника, определяющий его положение в волноводе для разных вычислительных экспериментов (см. рисунок 3); т - ширина во временной и частотной области; Az, Ax - шаги дискретизации, равномерные приращения по 2-D пространству; NC = 80 - количество узлов расчетной сетки на одну длину волны.
Так как для нашей задачи, в (1) мы упростили J и M независимые источники энергии как векторы плотности, то источники будут задаваться определёнными значениями полей напряженности в конкретные моменты времени, в соответствующих под областях тех или иных компонент поля. Например: Ex(i',k') = Spos(n), где i', k' - индексы определяющие области возбуждения (источника), где (i',k') некое подмножество множества (i,k) массива индексов определяющих Ex в данный момент времени.
Характер и геометрия неоднородности n_pml_xn_s на рисунке 3, т.е. нижней границы волновода Земля-Ионосфера, будет определяться в виде: 1-го приближения как на рисунке 2, и 2-го приближения, рисунок 4.
О 500 1000 1500 2000 2400
riz
Рис. 4. Профиль проводимости cre для неоднородности 2-м приближении [Figure 4. Conductivity profile cre for the inhomogeneity of the 2nd approximation ]
Результаты моделирования распространения атмосферика
Предполагается, что анализируя параметры атмосферика можно установить наличие неоднородности на трассе волновода, с целью чего проведено 6 компьютерных симуляций с двумя типами неоднородностей. В данных численных экспериментах варьировались 3 положения источника излучения как на рисунке 3. Далее представлены стоп-кадры каждой из симуляций в один и тот же момент времени 2200 отсчётов, где на рисунках (5-7, верх) неоднородность 1-го приближения, согласно рисунку 2; А на рисунках (5-7, низ) неоднородность 2-го приближения, согласно рисунку 4.
Заключение
С помощью математического моделирования и компьютерных симуляций процесса взаимодействия ЭМ волны с неоднородностью нижней границы (земля)
Mode - ТЕу I Display - Ez | Source - <~Ex(l Д40) ,WL=40000[m;|, 7500[Hz]>
• 10s time Step = 2200, time = 0.0025927 Sec ■ 10"7
Mode - ТЕу I Display - Ez I Source - <~ExC1.140),WL=40000[ffl].7500[Hz]>
¡105 time step = 2200, time = 0.0025927 sec ■ 10"'
г [meters!
Рис. 5. Результат компьютерной симуляции ЭМ излучения с источником S3(t) [Figure 5. Result of computer simulation of EM radiation with source S3(t) ]
Mode - TEy I Display - Ez I Source - <~Ex(l,300},WL.40000[m].7500[Hz]>
I - Um« -if-fj - 220(3, tin'* 0.0025Я 2 7 >
; rrii&C&rtl
Рис. 6. Результат компьютерной симуляции ЭМ излучения с источником S2(t) [Figure 6. Result of computer simulation of EM radiation with source S2(t) ]
Mode - TEy I Display - Ez I Source - <~£x(l,460),WL-40000[m].7500[Hz]>
1Ü5 time step = 2200. tints-0.002592? sec l'J '
6 8 10 5 z Tmetersl x 10
Mode - TEy | Display - Ez I Source - <~Ex(l,460),WL=40000[m],7500[Hz]>
X105 time step = 2200, time = 0.0025927 sec x 101
z Tmetersl
Рис. 7. Результат компьютерной симуляции ЭМ излучения с источником Si(t) [Figure 7. Result of computer simulation of EM radiation with source Si(t) ]
волновода, показано, что действительно имеет место обратное рассеяние ЭМ волны на трассе волновода. Обратное рассеяние возникает в следствии отражения ЭМ волны, при её взаимодействии с неоднородностью проводимости нижней границы волновода. Для разных приближений неоднородностей проводимости наблюдаются волны обратного рассеяния с разной интенсивностью и разным временем прихода.
Действительно, математическое и компьютерное моделирование процесса может помочь в установлении связи между параметрами излучения и неоднородностей. Можно предположить что верно и обратное: наблюдение за параметрами ЭМ сигнала определяющими атмосферик, в том числе наблюдаемой иногда волной обратного рассеяния, позволяет установить наличие неоднородности проводимости на трассе распространения.
Исходя из того что, такое локальное изменение проводимости может быть вызвано изменением химического состава подземных подповерхностных вод в следствии геологи-ческой активности предшествующей землетрясениям. Становится возможным, в будущем проверить гипотезу о том, что есть связь между параметрами атмосферика и некоторыми гидрогеохимическими предвестниками землетрясений.
Аббревиатуры
WWLLN The World Wide Lightning Location Network
VLF Very Low Frequency
EM Electromagnetic
ЭМ Электромагнитный
FDTD Finite-Difference Time-domain
TMy Transverse Magnetic to y
TEy Transverse Electric to y
PDE Partial differential equation
PEC Perfect electric conductor
ABC Absorbing boundary conditions
PML Perfectly matched layer
Список литературы
1. Koronczay, D., Lichtenberger, J.,Clilverd, M. A., Rodger, C. J., Lotz, S.I., Sannikov, D. V., et al. The source regions of whistlers., Journal of Geophysical Research: SpacePhysics, 2019. vol. 124, pp. 5082-5096 DOI: 10.1029/2019JA026559.
2. Lichtenberger J., Ferencz C., Bodnar L., Hamar D., Steinbach P.Automatic whistler detector and analyzer system: Automatic whistler detector, J. Geophys. Res., 2008. vol.113, no.A12 DOI: 10.1029/2008JA013467.
3. Storey L. R. O. An investigation of whistling atmospherics, Philosophical Transactions of the Royal Society of London. Series A, Mathematical and Physical Sciences, 1953. vol. 246, no. 908, pp. 113141 DOI: 10.1098/rsta.1953.0011.
4. Альперт Я. Л., Гусева Э.Г., Флигель Д. С. Распространение низкочастотных электромагнитных волн в волноводе Земля - ионосфера. Москва: Наука, 1967.123 с.
5. Budden K. G., Eve M. Degenerate modes in the Earth-ionosphere waveguide, Proceedings of The Royal Society A: Mathematical, Physical and Engineering Sciences, 1975. vol.342, no. 1629, pp. 175-190 DOI: 10.1098/rspa.1975.0019.
6. Milton K., Schwinger J. Electromagnetic Radiation: Variational Methods, Waveguides and Accelerators. New York: Springer, 2006. 357 pp. ISBN 978-3-540-29304-0 DOI: 10.1007/3-540-29306-X.
7. Sarkar T.K., Mailloux R. J., Oliner A. A., Palma M.S., Sengupta D. L. The Evolution of Electromagnetic Waveguides: From Hollow Metallic Guides to Microwave Integrated Circuits / History of Wireless, chapter 16. New-York, John Wiley & Sons, 2006, pp. 543-566 DOI: 10.1002/0471783021.ch16.
8. Spies K. P., Wait J. R. Mode calculations for VLF propagation in the earth-ionosphere waveguide. Washington, USA: U.S. Dept. of Commerce, 1961.116 pp.
9. Wait J.R. Electromagnetic Waves in Stratified Media. Oxford, UK: Oxford University Press, 1970. 384 pp. ISBN 9780198592235.
10. Helliwell R. A. Whistlers and Related Ionospheric Phenomena. New York: Dover Publications, 2006. 368 pp. ISBN 978-0486445724.
11. Helliwell R. A., Pytte A. Whistlers and Related Ionospheric Phenomena, American Journal of Physics, 1966. vol.34, no. 1, pp. 81-81 DOI: 10.1119/1.1972800.
12. Field Jr E. C., Gayer S.J. Dambrosio B. P. ELF propagation in the presence of nonstratified ionospheric disturbances, Final Report, May 1979 - Feb. 1980 Pacific-Sierra Research Corp., Santa Monica, CA, 1980. vol. 1, no. 1.
13. Lysak R. L., Yoshikawa A. Resonant Cavities and Waveguides in the Ionosphere and Atmosphere / : Magnetospheric ULF Waves: Synthesis and New Directions, American Geophysical Union, 2006, pp. 289-306 DOI: 10.1029/169GM19.
14. Kryukovsky A., Kutuza B., Stasevich V., Rastyagaev D. Ionospheric Inhomogeneities and Their Influences on the Earth's Remote Sensing from Space, Remote Sensing, 2022. vol.14, no. 21:5469, pp. 1-20 DOI: 10.3390/rs14215469.
15. Копылова Г. Н., Гусева Н.В., Копылова Ю.Г., Болдина С. В. Химический состав подземных вод режимных водопроявлений Петропавловского геодинамического полигона, Камчатка: типизация и эффекты сильных землетрясений, Journal of Volcanology and Seismology, 2018. №4, С. 43-62 DOI: 10.1134/S0203030618040041.
16. Chiodini G., Cardellini C., Di Luccio F., Selva J., Frondini F., Caliro, S., Rosiello A., Beddini G., Ventura G. Correlation between tectonic CO2 Earth degassing and seismicity is revealed by a 10-year record in the Apennines, Italy, Science Advances, 2020. vol. 6, no. 35, pp. e623 DOI: 10.1126/sciadv.abc2938.
17. Maxwell J. C. A dynamical theory of the electromagnetic field, Philosophical Transactions of the Royal Society of London, 1865. vol.155, pp. 459-512 DOI: 10.1098/rstl.1865.0008.
18. Jackson J.D. Classical Electrodynamics. 3rd Edition. New-York, USA: Wiley, 1998.832 с. ISBN 978-0471309321.
19. Yee K. Numerical solution of initial boundary value problems involving maxwell's equations in isotropic media, IEEE Transactions on Antennas and Propagation, 1966. vol. 14, no. 3, pp. 302307 DOI: 10.1109/TAP.1966.1138693.
20. Wang Z., Zhou C., Zhao S., Xu X., Liu M., Liu Y., Liao L., Shen X. Numerical Study of Global ELF Electromagnetic Wave Propagation with Respect to Lithosphere-Atmosphere-Ionosphere Coupling, Remote Sensing, 2021. vol.13, no. 20:4107, pp. 1-27 DOI: 10.3390/rs13204107.
21. Taflove A., Hagness S.C. Computational Electrodynamics: The Finite-difference Time-domain Method. Massachusetts, USA: Artech House, 2005. 1038 pp. ISBN 9781580538329.
22. Elsherbeni A.Z., Demir V. The finite-difference time-domain method for electromagnetics with MATLAB simulations. Raleigh, USA: SciTech Publishing, 2015. 560 pp. ISBN 978-1-61353-175-4.
23. Tverdyi D.A., Malkin E. I., Parovik R. I. Mathematical modeling of the propagation of a plane electromagnetic wave in a strip waveguide with inhomogeneous boundary conductivity, Bulletin KRASEC. Physical and Mathematical Sciences, 2022. Т. 41, №4, С. 66-88 DOI: 10.26117/20796641-2022-41-4-66-88.
24. Lindell I. V., Sihvola A. H. Perfect Electromagnetic Conductor, Journal of Electromagnetic Waves and Applications, 2005. vol.19, no. 7, pp. 861-869 DOI: 10.1163/156939305775468741.
25. Berenger J. P. A perfectly matched layer for the absorption of electromagnetic waves, Journal of Computational Physics, 1994. vol.114, no. 2, pp. 185-200 DOI: 10.1006/jcph.l994.1159.
26. Nickelson L. Electromagnetic Theory and Plasmonics for Engineers. Singapore: Springer, 2018.749 pp. ISBN 9789811323522.
27. Berenger J. P. Perfectly matched layer for the FDTD solution of wave-structure interaction problems, IEEE Transactions on Antennas and Propagation, 1996. vol. 44, no. 1, pp. 110-117 DOI: 10.1109/8.477535.
28. Andrew W. V., Balanis C.A. Tirkas P. A. Comparison of the Berenger perfectly matched layer and the Lindman higher-order ABC's for the FDTD method, IEEE Microwave and Guided Wave Letters, 1995. vol.5, no. 6, pp. 192-194 DOI: 10.1109/75.386128.
29. Veihl J. C., Mittra R. Efficient implementation of Berenger's perfectly matched layer (PML) for finite-difference time-domain mesh truncation, IEEE Transactions on Antennas and Propagation, 1996. vol.6, no. 2, pp. 94-96 DOI: 10.1109/75.482000.
Информация об авторах
Твёрдый Дмитрий АлександровичА - кандидат физико-математических наук, научный сотрудник лаборатории электромагнитных излучений института космофизических исследований и распространения радиоволн ДВО РАН, Паратунка, Россия, © СЖСГО 0000-0001-6983-5258.
Малкин Евгений ИльичА - научный сотрудник лаборатории электромагнитных излучений института космофизических исследований и распространения радиоволн ДВО РАН, Паратунка, Россия, © СЖСГО 0000-0001-8037-1335.
References
[1] Koronczay D., Lichtenberger J., Clilverd M. A., Rodger C. J., Lotz S.I., Sannikov D. V., et al. The source regions of whistlers, Journal of Geophysical Research: SpacePhysics, 2019, vol. 124, pp. 5082-5096. DOI: 10.1029/2019JA026559.
[2] Lichtenberger J., Ferencz C., Bodnar L., Hamar D., Steinbach P. Automatic whistler detector and analyzer system: Automatic whistler detector, J. Geophys. Res., 2008, vol. 113, no. A12. DOI: 10.1029/2008JA01346.
[3] Storey L.R. O. An investigation of whistling atmospherics, Philosophical Transactions of the Royal Society of London. Series A, Mathematical and Physical Sciences, 1953, vol. 246, no. 908, pp. 113-141. DOI: 10.1098/rsta.1953.0011.
[4] Alpert Y. L., Guseva E. G., Fligel D. S. Rasprostranenie nizkochastotnyh elektromagnitnyh voln v volnovode Zemlya - ionosfera [Propagation of low-frequency electromagnetic waves in the waveguide Earth - ionosphere]. Moscow, Science, 1967, 123 p.(In Russian).
[5] Budden K.G., Eve M. Degenerate modes in the Earth-ionosphere waveguide, Proceedings of The Royal Society A: Mathematical, Physical and Engineering Sciences, 1975, vol. 342, no. 1629, pp. 175-190. DOI: 10.1098/rspa.1975.0019.
[6] Milton K., Schwinger J. Electromagnetic Radiation: Variational Methods, Waveguides and Accelerators. New York, Springer, 2006, 357 pp. DOI: 10.1007/3-540-29306-X.
[7] Sarkar T.K., Mailloux R. J., Oliner A.A., Palma M.S., Sengupta D.L. The Evolution of Electromagnetic Waveguides: From Hollow Metallic Guides to Microwave Integrated Circuits, inbook: History of Wireless, chapter 16, 2006, pp. 543-566. DOI: 10.1002/0471783021.ch16.
[8] Spies K. P., Wait J. R. Mode calculations for VLF propagation in the earth-ionosphere waveguide. Washington, U.S. Dept. of Commerce, 1961, 116 p.
[9] Wait J. R. Electromagnetic Waves in Stratified Media. Oxford, UK, Oxford University Press, 1970, 384 p.
[10] Helliwell R. A. Whistlers and Related Ionospheric Phenomena. New York, Dover Publications, 2006, 368 p.
[11] Helliwell R. A., Pytte A. Whistlers and Related Ionospheric Phenomena, American Journal of Physics, 1966, vol. 34, no. 1, pp. 81-81. DOI: 10.1119/1.1972800.
[12] Field Jr E. C., Gayer S. J. Dambrosio B. P. ELF propagation in the presence of nonstratified ionospheric disturbances, Final Report, May 1979 - Feb. 1980 Pacific-Sierra Research Corp., Santa Monica, CA., 1980, vol. 1, no. 1.
[13] Lysak R. L., Yoshikawa A. Resonant Cavities and Waveguides in the Ionosphere and Atmosphere, In book: Magnetospheric ULF Waves, American Geophysical Union, 2006, pp. 289-306, DOI: 10.1029/169GM19.
[14] Kryukovsky A., Kutuza B., Stasevich V., Rastyagaev D. Ionospheric Inhomogeneities and Their Influences on the Earth's Remote Sensing from Space, Remote Sensing, 2022, vol. 14, no. 21:5469, pp. 1-20. DOI: 10.3390/rs14215469.
[15] Kopylova G. N., Guseva N.V., Kopylova Y. G., Boldina S.V. Himicheskij sostav podzemnyh vod rezhimnyh vodoproyavlenij Petropavlovskogo geodinamicheskogo poligona, Kamchatka: tipizaciya i effekty sil'nyh zemletryasenij [Chemical composition of groundwater in regime water manifestations of the Petropavlovsk geodynamic test site, Kamchatka: typification and effects of strong earthquakes], Journal of Volcanology and Seismology, 2018, no. 4, pp. 43-62. DOI: 10.1134/S0203030618040041.
[16] Chiodini G., Cardellini C., Di Luccio F., Selva J., Frondini F., Caliro, S., Rosiello A., Beddini G., Ventura G. Correlation between tectonic CO2 Earth degassing and seismicity
is revealed by a 10-year record in the Apennines, Italy, Science Advances, 2020, no. 35, pp. e623. DOI: 10.1126/sciadv.abc2938.
[17] Maxwell J. C. A dynamical theory of the electromagnetic field, Philosophical Transactions of the Royal Society of London, 1865, vol. 155, pp. 459-512. DOI: 10.1098/rstl.1865.0008.
[18] Jackson J. D. Classical Electrodynamics. 3rd Edition. New-York, USA, Wiley, 1998, 832 p.
[19] Yee K. Numerical solution of initial boundary value problems involving maxwell's equations in isotropic media, IEEE Transactions on Antennas and Propagation, 1966, vol. 14, no. 3, pp. 302-307. DOI: 10.1109/TAP.1966.1138693.
[20] Wang Z., Zhou C., Zhao S., Xu X., Liu M., Liu Y., Liao L., Shen X. Numerical Study of Global ELF Electromagnetic Wave Propagation with Respect to Lithosphere-Atmosphere-Ionosphere Coupling, Remote Sensing, 2021, vol. 13, no. 20:4107, pp. 1-27. DOI: 10.3390/rs13204107.
[21] Taflove A., Hagness S. C. Computational Electrodynamics: The Finite-difference Timedomain Method. Massachusetts, USA, Artech House, 2005, 1038 p.
[22] Elsherbeni A. Z., Demir V. The finite-difference time-domain method for electromagnetics with MATLAB simulations. Raleigh, USA, SciTech Publishing, 2015, 560 p.
[23] Tverdyi D.A., Malkin E. I., Parovik R.I. Mathematical modeling of the propagation of a plane electromagnetic wave in a strip waveguide with inhomogeneous boundary conductivity, Bulletin KRASEC. Phys. and Math. Sci. 2022, vol. 41, no. 4, pp. 66-88. DOI: 10.26117/2079-6641-2022-41-4-66-88.
[24] Lindell I.V., Sihvola A. H. Perfect Electromagnetic Conductor, Journal of Electromagnetic Waves and Applications, 2005, vol. 19, no. 7, pp. 861-869. DOI: 10.1163/156939305775468741.
[25] Berenger J. P. A perfectly matched layer for the absorption of electromagnetic waves, Journal of Computational Physics, 1994, vol. 114, no. 2, pp. 185-200. DOI: 10.1006/jcph.1994.1159.
[26] Nickelson L. Electromagnetic Theory and Plasmonics for Engineers. Singapore, Springer, 2018, 749 p.
[27] Berenger J. P. Perfectly matched layer for the FDTD solution of wave-structure interaction problems, IEEE Transactions on Antennas and Propagation, 1996, vol. 44, no. 1, pp. 110117. DOI: 10.1109/8.477535.
[28] Andrew W. V., Balanis C.A. Tirkas P. A. Comparison of the Berenger perfectly matched layer and the Lindman higher-order ABC's for the FDTD method, IEEE Microwave and Guided Wave Letters, 1995, vol. 5, no. 6, pp. 192-194. DOI: 10.1109/75.386128.
[29] Veihl J. C., Mittra R. Efficient implementation of Berenger's perfectly matched layer (PML) for finite-difference time-domain mesh truncation, IEEE Transactions on Antennas and Propagation, 1996, vol. 6, no. 2, pp. 94-96. DOI: 10.1109/75.482000.
Information about authors
Tverdyi Dmitrii AlexsandrovichA - Ph. D. (Phys. & Math.), Researcher, laboratory of electromagnetic propogation Institute of Cosmophysical Research and Radio Wave Propagation FEB RAS, Paratunka, Russia, <E ORCID 0000-0001-6983-5258.
Malkin Evgeniy IlichA - Researcher, Laboratory of Electromagnetic Radiation Institute of Cosmophysical Research and Radio Wave Propagation FEB RAS, Paratunka, Russia, ©ORCID 0000-0001-8037-1335.
