Численное моделирование распространения электромагнитных волн СДВ-диапазона в волноводном канале Земля-ионосфера в ближней области источника Текст научной статьи по специальности «Физика»

В. М. Ромаданов, Л. Н. Лутченко

ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ РАСПРОСТРАНЕНИЯ ЭЛЕКТРОМАГНИТНЫХ ВОЛН СДВ-ДИАПАЗОНА В ВОЛНОВОДНОМ КАНАЛЕ ЗЕМЛЯ—ИОНОСФЕРА В БЛИЖНЕЙ ОБЛАСТИ ИСТОЧНИКА

Введение. Современные пакеты для численного моделирования, основанные на методе конечных элементов, использующие алгоритмы автоматической генерации сетки, становятся неэффективными с ростом количества узлов сетки. Для работы универсальных алгоритмов, которые строят сетку и составляют матрицы, применяемые в методе конечных элементов, необходимо больше памяти, чем для численного решения самой системы уравнений, если элементы матрицы уже найдены. Но если для решения задачи не требуются ресурсоёмкие алгоритмы, то решение системы уравнений может быть существенно оптимизировано за счёт построения регулярной сетки. При этом сетка может быть неравномерной, что не приводит к существенному усложнению алгоритмов. Численные методы, применяющиеся для решения уравнений электродинамики, описаны в [1]. В [2] предложен приём, позволяющий моделировать распространение электромагнитных волн в неограниченном пространстве, который не приводит к существенному увеличению расхода времени и памяти, необходимых для численного решения уравнений.

Постановка задачи. Решается система уравнений Максвелла в цилиндрической системе координат (r, ф^). Зависимость от времени считается гармонической: exp(-irot). Система уравнений Максвелла для комплексных амплитуд полей:

rot(E) = iro^0H, (1)

rot (iH) = -Ш£0£.Ё, (2)

где ro = 2nf; f = 11, 9 кГц — частота; c — скорость света в вакууме; е — относительная комплексная диэлектрическая проницаемость ионосферы, которая может зависеть от координат r и z. Считаем, что в волноводе возбуждается только одна TM-поляризация, тогда электрическое и магнитное поля имеют компоненты Ё = (Er, 0,Ez) и Н = (0, Нф, 0) соответственно. На поверхности Земли (z = 0) ставятся импедансные граничные условия dH

—^ + *А-6Яф|г=п = 0, (3)

uz z=0

где k = ro/c, 5 — приведённый поверхностный импеданс подстилающей поверхности. Сверху на высоте zPML, а также на некотором удалении грмь от источника располагаются слои со специально подобранными свойствами среды так, чтобы в них происходило поглощение электромагнитных волн. На внешних границах этих слоёв ставится условие, соответствующее идеально проводящей стенке, как показано на рис. 1:

Ez\r=r = 0 Er\z=z = ° (4)

4 —' max 'z — zmax 4 '

Уравнения (1), (2) однородные и справедливы в области, свободной от источников, поэтому внутренняя граница области, где ищется численное решение, выбирается на

© В. М. Ромаданов, Л. Н. Лутченко, 2011

Рис. 1. Область, в которой ищется решение и слои для поглощения электромагнитных волн, уходящих за границы этой области

E = 0

r

PML PML

r 1 r 1 1 r s

z = 4V/)1'2 z • = z^/)1'2

PML

e'Jrz) 1 r

р = 1 z = z

/

/

/

/

/

/

/

/

/

/

/

'/ °

/ II

/ г ^

/

/

/

/

/

/

/

/

/

/

/

/

r r

max

max

Z

PML

0

PML

расстоянии rm;n от оси, на которой расположен источник. На этой границе задаётся комплексная амплитуда поля.

Расстояние rmin меньше длины волны в вакууме, но больше длины волны в земле. Внешняя граница выбирается на удалении от источника rmax, много большем длины волны, где волна распространяется в направлении от источника. На внешней стенке используются неотражающие граничные условия [2, 3]. При поиске аналитического решения в неограниченной среде используют специальные условия на бесконечности, чтобы получить в решении волны, распространяющиеся от источника. Для численного моделирования необходимо ограничить область, в которой ищется решение, а на внешней стенке поставить условия, при которых электромагнитные волны проходят через эту границу без отражения так, как будто они распространяются в неограниченном пространстве. Один из вариантов такого граничного условия называется Perfectly Matched Layer (PML).

Относительная комплексная диэлектрическая проницаемость ионосферы е в (2) в случае регулярного волновода зависит только от вертикальной координаты z. Зависимость выражается формулой [4]

e(z) = 1 -

4jt2 80,6 • l06Ne(z) ю(ю + iv(z))

(5)

где функция Ые(г) — описывает концентрацию электронов на кубический сантиметр (рис. 2); у(г) — частота соударения электронов с нейтральными частицами,

у(г) = 45,45 • 105 ехр(0,148(70 - г)).

Неотражающие граничные условия.

Для численного решения уравнений Максвелла необходимо построить сетку, которая всегда имеет ограниченное количество узлов. Кроме того, с целью сокращения расхода памяти компьютера и времени на

Рис. 2. Плотность электронной концентрации

решение уравнений, размеры области решения необходимо выбирать как можно меньше. Таким образом, при численном решении задач, в которых волны могут распространяться в неограниченном пространстве, требуется ввести искусственную границу, отсекающую внешнюю часть задачи, но при этом обеспечивающую прохождение электромагнитных волн так, как будто они распространяются в неограниченном пространстве. Одним из вариантов такого типа граничных условий является PML-ABC (Perfectly Matched Layer Adsorbed Boundary Condition), предложенный Беренгером [2, 3]. Неотражающая граница состоит из идеально проводящей стенки, перед которой расположен слой среды со специальными свойствами. В этой среде происходит затухание волн, но при этом нет отражения и преломления. Электромагнитные волны PML описываются модифицированными уравнениями Максвелла, в которые добавлены дополнительные проводимости так, чтобы оба уравнения выглядели похожим образом.

Для выбора параметров, описывающих поглощающий слой, рассмотрим более простую задачу, в которой в цилиндрической системе координат волна падает на расположенную на расстоянии го от источника вертикальную границу между обычной средой с комплексной относительной диэлектрической проницаемостью £1 и специальной средой с проницаемостью £2 и дополнительными проводимостями 02 и о2.

Распишем уравнения (1) и (2) по компонентам для первой среды

дЕг дЕ~~ ■ ТТ

-qZ- - = *«¥о#ф, (6)

дНф

= иое0е1Ег, (7)

1 д

(гЯф) = -iw£0£iEz, (8)

и для второй среды по аналогии с тем, как это сделано в [2, 3],

дЕг

~qZ~ (9)

дЕ

дг

дНф

— $Ш^0.Нфг , (10)

^ = *Ш£0£1 Ег, (И)

1 д

— — (?\Нф) = —iш£■Q£■lEz + ооЕг = —*ш[1о вг-Бг, (12)

где Яф Нфг + . Да.лее

= 1 + -^-, 4 = 1 + —. (13)

Ш8082 Ш^о

Из уравнений (7) и (8) выражаются компоненты электрического поля

„ _ 1 дЯФ

r *Ю£о£1 dz

Ей = -^—~{гЩ).

iW£o£l г дг

которые подставляются в (6). В результате получается уравнение для магнитного поля

Общее решение уравнения (15) может быть выражено в виде суммы полей падающей и отражённой волн

жённую волну; Яо — произвольная постоянная, которая имеет смысл амплитуды магнитного поля падающей волны, компоненты волнового вектора подчиняются дисперсионному соотношению к2г + к^ — к2£1; К — коэффициент отражения. Соответственно из (16) и (14) получается формула для вертикальной компоненты электрического поля слева от границы

В области PML, расположенной справа от границы, из (9)—(12) выражаются компоненты электрического поля

Возникшие коэффициенты л/вгв* в дифференциальном операторе перед производными по г можно интерпретировать как коэффициенты перехода в другую систему координат [2]. Если выполнить замену переменных г' = у/вгв*г, уравнения (17) и (15) будут различаться только значением коэффициента в последнем слагаемом. Справа от границы решение уравнения должно быть в виде прошедшей волны

где для компонент выполняется дисперсионное соотношение к2,г + к|2 — к2 £2; Т — коэффициент прохождения.

На границе двух сред касательные компоненты Яф и Ег должны быть непрерывны, а также должно выполняться равенство кг1 — кг2. С этими условиями получаются уравнения

(15)

где к2 = ю2£о^о-

(16)

где нО^ (к1гг) и Я(1) (к1гг) — функции Ханкеля первого рода, описывающие падающую

(2) (2)

волну; Я0 (к1Гг) и Я} (к1Гг) — функции Ханкеля второго рода, описывающие отра-

и уравнение для магнитного поля

(17)

#!(1) (ku.ro) + нн[2) (к1гг0) = ТН^ (к2Гу/^*го) ,

из которых находятся К и Т:

в г в* г о) ГЯХ(1) (А;г1/«г«г?’о)

•гПз) ^ Я^1) (Ауг0)

--------------Л/ 5-Г 77\-------------------

Я} ; (кг^вгвр'о) Я} > (кт^вгвр'о)

Для упрощения (18) воспользуемся асимптотическим представлением функций Ханке-

В результате приближенное выражение для коэффициента отражения примет вид

Чтобы Н = О, должно выполняться условие вг = в* откуда следует равенство = —

£0£2 Ц0

в соответствии с введёнными определениями (13).

Выбрав £1 = £2, получаем в обеих средах одинаковое дисперсионное соотношение, следовательно, компонента ки = к2Г = кг прошедшей волны будет такой же, как и у падающей, потому что на границе компонента кг должна быть непрерывной. Из условия вг = в* следует, что л/вгв* = вг. Компонента магнитного поля справа от границы Нф ~ Я(1) (гкгвгг). Предположим, что волна распространяется в направлении вдоль оси г, тогда кг = Определим толщину слоя РМЬ, в котором происходит затуха-

ние, в е раз:

Параметр 02 отвечает за скорость затухания.

Для численного моделирования распространения электромагнитных волн в неограниченной среде используется слой конечной толщины с параметрами (13), поэтому из этого слоя в обычную среду выйдет волна, отражённая от внешней границы. Параметр

02 необходимо выбирать так, чтобы за счёт затухания эта отражённая волна была пренебрежимо малой амплитуды.

Численная схема. Решение описанной выше задачи выполняется методом конечных разностей. Причём используется неравномерная регулярная прямоугольная сетка. Выбор в пользу этого метода сделан из следующих соображений. Использование конечно-элементных алгоритмов, в которых происходит автоматическое создание сетки,

у цою/ с

— Іт [^] + 11е [у/го\ ) г ехр г

с ростом количества узлов сетки оказывается неэффективным. Процесс создания сетки и вычисления элементов матриц для метода конечных элементов требует заметно больше аппаратных ресурсов по сравнению с алгоритмами решения систем линейных алгебраических уравнений с разрежённой матрицей. Если использовать регулярную сетку и самостоятельно определить связь между положением узла в сетке и его порядковым номером, то можно существенно сократить расход памяти компьютера и время счёта.

Для решения уравнений в прямоугольной области используется прямоугольная сетка размером М узлов по горизонтали в направлении оси г и N узлов по вертикали в направлении оси г. Для численного решения уравнений в частных производных сеточным методом необходимо пронумеровать узлы сетки. В методе конечных разностей дифференциальное уравнение заменяется линейной системой алгебраических уравнений. Удобно ввести нумерацию узлов сетки так, чтобы номер узла соответствовал номеру элемента в векторе неизвестных в системе алгебраических уравнений. Нет необходимости нумеровать каждый узел сетки. Во-первых, узлы, расположенные по углам прямоугольной области, не используются вообще. Во-вторых, в узлах, расположенных в вертикальном сечении на границе г = гт;п, комплексная амплитуда поля Нф задана, поэтому их необходимо пронумеровать отдельно от остальных узлов, так как эти узлы относятся к вектору в правой части системы алгебраических уравнений. Пусть нумерация узлов начинается с нуля. Индексы т и п обозначают положение узла в сетке и соответствуют координатам г и г. Номер кт,п соответствует номеру элемента в векторе неизвестных. Связь кт,п с т и п определена формулой

^N (т — 1) + п, т<М — 1

кт,п — Л ,г/ ч , (19)

(т — 1) + п — 1, т = М — 1.

Количество элементов в векторе неизвестных системы алгебраических уравнений равно (М — 1^ — 2, а количество ненулевых элементов в разрежённой матрице находится по формуле 5(М — 3)^ — 2) +4^ — 2) + 2((М — 2)2 + N — 2). Обратный переход от номера кт,п к индексам т и п осуществляется по формулам

к

N

+ 1, (20)

|к — N (т — 1), т<М — 1

п = (21)

|^к — N (т — 1) + 1, т = М — 1.

Производные аппроксимировались конечными разностями на неравномерной сетке по формулам

■7^7 = (°1^£-1 + 02Ц; + °3^*+1) 1 (22)

д2и 1

0^2 = ~Ц) + ^2 и г + Ь^и^х) , (23)

где вместо и может быть либо Нф, либо £, вместо х — либо г, либо г, индекс г — либо

т, либо п. Коэффициенты вычисляются по координатам узлов

П = — [(*<_! — ж*)(ж*+! — ХгУ — (ж*-)-! — ж*)(ж*_1 — ХгУ] , (24)

О! = ^(ж*+1 - ж*)2, (25)

02 = \{хг+1 - ж*_1)(-ж*_1 + 2ж* - Ж*+1), (26)

03 = — Хг—\)", (27)

Ъ\ = —х^+1 + х*, (28)

Ь2 = х*+1 — х_1, (29)

Ъз = х_1 — хн. (30)

Из уравнений (1), (2) получается дифференциальное уравнение для комплексной амплитуды компоненты магнитного поля Нф

1 а / знл глн* _ 1 л 1 ^ _ 1 а» ан» , =

V У 5 г 2 е <9г г <9г г дг дг

г 9/’ у 9/’ у дг2 е дг г дг е дг дг ^

Это уравнение может быть записано в виде

д2Щ , д2Щ , дЩ , дЩ , ^ _п ,о1Х

01 “а?7" С2_аГ ф ““ ’ (31)

где коэффициенты

11 де 1 де 2 1 д е

С1 =-------ТТ-! с2 =-7Г-1 С3 = АГе-------—.

г е дг е дг ге дг

Заменив производные в (31) конечными разностями в соответствии с (22) и (23), получаем систему линейных алгебраических уравнений

^1Ят_1,п + ^2Ят+1,п + ^3 Ят,п_1 + ^4 Ят,п+1 + ^5Ят,п = 0 (32)

с коэффициентами, которые также зависят от пары индексов т, п:

1

Й2 = ^; йз = ж

4 = ^

Ъу 1 + С1ЯГ1).

Ъг3 + с1агз),

Ъг1 + С2021),

Ъг3 + с2ог3):

Ьг2 + С10г2) + у— (Ьг2 + с2 0-2) + Сз-

^ г

г -‘-''г

Обозначение означает, что в (24) вместо х подставлено г, а вместо г — т. Аналогично означает, что в (24) вместо х подставлено г, а вместо г — п. Точно такой же принцип используется при обозначениях оГ1, оГ2, огз, ог1, ог2, огз, ЪГ1, ЪГ2, Ъгз, Ъг1, Ъг2, Ъгз для

коэффициентов (25)-(30).

Из граничного условия (3), если заменить производные конечными разностями, получаются разностные уравнения для значений Яф в узлах, расположенных на нижней границе,

~Ят,1 + \ikbm----;----— ) Нт,0 = 0. (33)

г1 — г0 V г1 — го

Граничные условия (4) эквивалентны

г дг^’Н^

дНф

дг

Заменив производные конечными разностями, получаем

1

Гм-1 - Гм-2

~НМ-2,п +

1

+

1

Нм

М-1,-

0.

Гм-1 — Гм-2 Гм-1,

— Нт,Ы-2 + Нт,Ы-1 = °

(34)

(35)

При составлении матрицы системы алгебраических уравнений последовательно перебираются все узлы сетки, в которых значение Нф неизвестно. Это все узлы, кроме угловых и узлов в начальном вертикальном сечении. Для каждого узла по формуле (19) определяется номер соответствующего ему элемента в векторе неизвестных. Этот же номер определяет и номер уравнения в системе. Соответственно коэффициент в (32) — это диагональный элемент. В (33) таким элементом является коэффициент перед Нт,о, в (34) — коэффициент перед Нм-1,п, в (35) — коэффициент перед Нт,н-1.

Далее для соседних узлов по формуле (19) определяется номер элемента в строке матрицы. Эти элементы соответствуют коэффициентам сІ1, сІ2, ^з, ^4 в (32). Если в процессе составления матрицы в (32) попадается узел с номером т = 0, то это слагаемое переносится в правую часть системы уравнений. Причём номер элемента вектора правой части совпадает с номером строки матрицы.

Чтобы определить, какому узлу сетки соответствует элемент с номером к вектора неизвестных в составленной системе уравнений, используются формулы (20) и (21).

Результаты моделирования. С помощью разработанной программы возможен расчёт комплексной амплитуды поля в волноводе, образованном поверхностью Земли и ионосферой, и над открытой поверхностью Земли. Выполнив расчёт комплексной амплитуды магнитного поля над открытой идеально проводящей поверхностью и либо над поверхностью с ненулевым импедансом, либо над поверхностью с ненулевым импедансом в волноводе, можно определить функцию ослабления (рис. 3).

По результатам численного моделирования также можно определить среднее за период значение вектора Умова—Пойнтинга

Мае

ЕхН *

Выражения для компонент Б после подстановки Е из (2) имеют вид

Бг = — 11е 2

1

*Ю£0Є 1

ІЮЄоЄ

Я* ф

Н

дЩ

ф дг

Заменив производные конечными разностями в соответствии с (22) и (24)—(27), получим формулы для определения значений Бг и в узлах сетки.

В результатах численного моделирования при всех рассмотренных импедансах было обнаружено, что амплитуда поля Нф имеет локальный минимум у поверхности Земли на расстоянии между 250 и 300 км от источника, а на расстоянии около 250 км от источника на высоте 30 км — локальный максимум.

г = г

z=z

\Щ arg(W),рад

Рис. 3. Зависимость функции ослабления от высоты в сечениях, расположенных до и после локального минимума

Интересно направление распространения энергии в этой области радиотрассы в окрестностях локальных минимумов и максимумов амплитуды поля. По результатам численного моделирования были найдены зависимости среднего значения за период компонент вектора Умова—Пойнтинга от координат в плоскости радиотрассы. На рис. 4 представлена зависимость вертикальной компоненты вектора Умова—Пойнтинга от высоты в двух вертикальных сечениях радиотрассы. Сечения расположены на расстояниях 250 и 300 км от источника, то есть по обе стороны от локального минимума амплитуды Нф. Из графиков видно, что на расстоянии 250 км от источника, на высотах между 10 и 20 км вертикальная компонента имеет экстремум и положительна, что соответствует направлению распространения энергии кверху. А на расстоянии 300 км от источника на тех же высотах компонента также имеет экстремум, но принимает отрицательные значения, что соответствует направлению распространения энергии книзу. Зависимость у самой поверхности Земли отрицательна, что выделено на рисунке. Это объясняется затуханием амплитуды Нф за счёт положительной вещественной части импеданса Земли.

На рис. 5 и 6 приведены зависимости функции ослабления от высоты в сечении, расположенном на расстоянии 1000 км от источника. Зависимости получены для сильно индуктивного импеданса, слабо индуктивного и ёмкостного. Качественное различие между зависимостями объясняется наличием полного внутреннего отражения волн СДВ диапазона от нижних слоёв ионосферы.

Вт/м

Рис. 4. Зависимость Яг от высоты в сечениях, расположенных до и после локального минимума

В [5] приводится аналитическое решение для функции ослабления поля излучения над открытой плоской поверхностью с поверхностным импедансом 5:

вш\|/соя2 \|/ 1 (1 -Ь2)Ь _ / г>л

^изл. — <£ , ... Н~~ ъ I ...

5 + Бт у 5 + Бт у

(ей) = 1 + 2%/вДехр( — аК) J ехр(г2)с1г,

■\Zs~R

в = гк(1 + 6эту — \/1 — Ь2 ссеу),

где у — угол скольжения; К — расстояние от источника до точки наблюдения. Если точка наблюдения расположена на поверхности Земли, то у = 0, и функция ослабления

^изл. = (1 - 52)ш(зЕ),

(36)

ш

в = ik( 1 — \/1 — Ь2).

ДОЯ (,2 ^ ДОЛ ‘2

ис

\Щ а^(^), рад

5. Зависимость функции ослабления от высоты в волноводе в сечении на расстоянии 1000 км от источника при разных импедансах Земли

\^| а^(^), рад

Рис. 6. Зависимость функции ослабления от высоты над открытой поверхностью в сечении на расстоянии 1000 км от источника при разных импедансах Земли

Результаты численного моделирования и аналитического решения (36) для функции ослабления на поверхности Земли на расстоянии 1000 км от источника при различных импедансах приведены в таблице.

Значения №изл. на поверхности Земли на расстоянии 1000 км от источника

при различных импедансах

5 Аналитическое Численное

0,0062 - 0,0094г 1,0393 ехр(0,2210г) 1,0388 ехр(0,2209г)

0,0062 - 0,0034г 0,9594 ехр(0,1329г) 0,9592 ехр(0,1327г)

0,0062 + 0,0034г 0,8762 ехр(0,0377г) 0,8764 ехр(0,0375г)

Выводы. Написана программа для расчёта комплексной амплитуды магнитного поля в волноводе и над открытой поверхностью Земли. Было проведено сравнение результатов численного решения с аналитическим решением для функции ослабления. Получено совпадение трёх знаков после запятой для модуля и для аргумента функции ослабления, при этом использовался шаг сетки 1 км в горизонтальном и вертикальном направлениях.

С помощью полученных результатов вычислено среднее значение за период вектора Умова—Пойнтинга и продемонстрировано, что перенос энергии в волноводе осуществляется в основном через те участки пространства, в которых амплитуда поля имеет максимумы в результате интерференции. Показано, что у поверхности Земли вертикальная составляющая вектора Умова—Пойнтинга направлена вниз.

Программа написана на C+—+ и оптимизирована по времени счёта и расходу памяти с учётом того, что для численного решения используется регулярная неравномерная сетка. Это позволило решить задачу, которую не удаётся решить на компьютере с теми же характеристиками с помощью известных коммерческих программ, использующих метод конечных элементов и автоматическую генерацию сетки.

Литература

1. Van Rienen U. Numerical Methods in Computational Electrodynamics. Berlin, 2001. 375 p.

2. Berenger J.-P. Perfectly Matched Layer (PML) for Computational Electromagnetics // Synthesis Lectures On Computational Electromagnetics. N 8. Morgan & Claypool Publisher. 2007.

3. Berenger J.-P. A Perfectly Matched Layer for The Absorption of Electromagnetic Waves // J. Comp. Phys. 1994. Vol. 114. P. 185-200.

4. Макаров Г. И., Новиков В. В., Рыбачек С. Т. Распространение электромагнитных волн над земной поверхностью. М., 1991. 196 с.

5. Новиков В. В. Распространение радиоволн над слоистой трассой // Проблемы дифракции и распространения волн. Вып. 1. Л., 1962. С. 116-132.

Статья поступила в редакцию 28 декабря 2010 г.