Численная схема решения уравнений Максвелла на примере гармонического электрического диполя Текст научной статьи по специальности «Физика»

ВЫЧИСЛИТЕЛЬНЫЕ МЕТОДЫ И ИНФОРМАЦИОННЫЕ ТЕХНОЛОГИИ В

ГЕОФИЗИКЕ

УДК 519.626

О. И. Ахметов, И. В. Мингалев, О. В. Мингалев, Ю. В. Федоренко

ЧИСЛЕННАЯ СХЕМА РЕШЕНИЯ УРАВНЕНИЙ МАКСВЕЛЛА НА ПРИМЕРЕ ГАРМОНИЧЕСКОГО ЭЛЕКТРИЧЕСКОГО ДИПОЛЯ

Аннотация

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

Ключевые слова:

распространение электромагнитных волн в простых средах, численное моделирование, гармонический диполь Герца.

O. I. Akhmetov, I. V. Mingalev, O. V. Mingalev, Yu. V. Fedorenko

DIFFERENCE SCHEME FOR THE NUMERICAL SOLUTION

OF MAXWELL’S EQUATIONS ON THE EXAMPLE OF HARMONIC ELECTRICAL

DIPOLE

Abstract

Test results of new difference scheme for the numerical solution of Maxwell's equations on spatial two dimensional analytical solution are presented for point source of field harmonic electrical Hertz dipole in infinite homogeneous isotropic medium, which is dielectric or conductor. The numerical solution is in good according to analytical dipole field in middle field and part of far field.

Keywords:

electromagnetic waves propagation in an ordinary medium, numerical modeling, harmonic Hertz dipole.

Введение

Начиная с 1994 г. за рубежом опубликовано большое число работ (см. обзоры [1, 2]), в которых для исследования распространения

электромагнитных сигналов в нижнем диапазоне частот используются различные численные модели, основанные на методе конечных разностей во временной области. В нашей работе [3] предложена численная модель с учетом тензорной проводимости ионосферы, основанная на численном решении системы уравнений Максвелла с законом Ома при помощи новой явной двухслойной по времени разностной схемы. В этой схеме используется

120

новый подход к аппроксимации по времени, который основан на представлении уравнений Максвелла в интегральной по времени форме.

Расчеты и результаты

В работе используются сферические координаты (r, а, в) с центром в источнике и осью вдоль единичного вектора направления момента диполя v (r=|x|, в - широта). Пусть ю, I0 и l - частота диполя, амплитуда тока и расстояние в нем, P0 = I01 /(4п) - амплитуда момента диполя в системе СИ.

Тестовые расчеты выполняются следующим образом. Для имитации источника поле от аналитического решения при t > 0 задается на сетке в шаре радиуса R0, причем плавно включается за время t0=10 0, где 0 =2п/ю - период диполя. То есть поле от аналитического решения умножается на множитель n(t /t0) - достаточно гладкую неубывающую функцию, такую, что p(s) = 0 при 5 < 0 и p(s) = 1 при 5 > 1. В сферическом слое {R0 < r < R1} поле рассчитывается путем численного решения уравнений Максвелла в однородной изотропной среде (е = s1 = Const, ц = s1 = Const, о = Const):

—=-rotE, —= rotH - j, divD = divB = 0, j = uE, D = ss0E, B = рр0Н. (1)

dt dt

Для имитации неограниченного пространства (то есть для подавления численного отражения от границы) на границе области моделирования {r = R1} используется вариант поглощающего слоя (UPML - uniaxial perfectly matched layer), когда поля рассчитываются на сетке в большем сферическом слое {R0 < r < R2, R2 > R1}. В поглощающем слое { R1 < r < R2} решаются уравнения (1) с возрастающими по r проницаемостями среды:

s=s(r) = s1 (l + yp((r-R)/(R2 -Ri))), P=p(r) = P1(l + yp((r-R1 )/(R2 -R1))),

где у > 0 - заданная постоянная, а функция p(s) описана выше.

В представленных расчетах использовалось у =19, т. е. проницаемости увеличиваются в поглощающем слое в у +1 = 20 раз. Это дает при r ^ R2

замедление скорости волны a (r) = c/^s(r)р(r) в 20 раз по сравнению

со скоростью в области моделирования a 1 = cj , а также уменьшает в разы

длину волны и приводит к уменьшению на каждом шаге величины распространяющихся в поглощающем слое полей. В результате в случае достаточно большого расстояния R2 — R1 поля практически обращаются в 0 при r ^ R2, когда из-за односторонних разностей при вычислении пространственных производных на внешней границе поглощающего слоя {r = R2} возникает численное отражение.

В случае диэлектрика (о = 0) дипольное аналитическое решение является точным решением уравнений Максвелла и в системе СИ имеет вид:

A( r,e,t ) = — e(at kr), H (r,P, t)= rot A (r,P, t)

E ( r,P,t ) =

rot H ( r,P, t)

l Ю SSo

(tot-kr)

l®SSo r"

(k 2r 2-1 -1

ikr)v-

= pvxx](i + ikr)el(*t-kr), (2)

r3

fcix(3 + 3ikr - k2r 2)j , (3)

121

где волновое число к = к0 (ю) = jс е . Размеры области моделирования

и поглощающего слоя были R1=4.3X, R2=4.8X, где X =2nlk - длина волны. Шаги пространственной сетки составляли Аг= X1 NX =Xl 160 и Ар = nl Np = nl360 рад = 0.5 град., т. е. число шагов было Nr=768, Np=360 и на поглощающий слой приходится 80 шагов по г. Шаг по времени т должен удовлетворять условию Куранта (сеточная скорость Vh больше скорости света в среде а1), которое в рассматриваемом случае удобно представить в виде:

minVh = min((R ДР/т); (Дг/т)) > a, = сЦs,^, . (4)

Рассматривались два варианта значения R0 и шага по времени т = 0lNe: R01=0.5X, Ti= 01720, Nei=720 и R02=1.5X, Т2= 01360, Ne2=360.

Рис.1. Зависимость для каждой компоненты поля рассчитанной по формуле (5)

усредненной установившейся относительной ошибки (§F (г, Р = 0})

от г в плоскости поперек диполя (Р=0): слева а - при R02=1.5X, т2= 0/360 и справа б - при R01=0.5X, т1= 0/720

На рисунке 1 для присутствующих компонент поля Ег, Ep, Ва представлены оценки зависимости от расстояния г до источника усредненной

установившейся относительной ошибки (§f (г,Р = 0}) численного решения

Fh (гр, р = 0, к)

по сравнению с аналитическим решением

Fa (гр, Р = 0, tk )

в плоскости поперек диполя (Р=0), (где F = Ег, Ер,Ва, Гр = рДг, tk = кт), которая рассчитывалась по формуле:

(

ку+К

\ If

(SF (гр, Р = 0}) = Z |Fh (гр , Р = 0, tk )-Fa (гр ,Р= 0, tk )|

, к=к +1

к,+К

Z Fa (

k=k +1

г, ,Р =

0,4

(5)

Усреднение (в установившемся режиме) проводилось по интервалу [Е1; T2], где T1=7000 и T2=19000, т. е. k1=700 N0, K=1200 N0. На рис.1а представлены результаты для R02=1.5X, т2= 61360, а на рис.1б при R01=0.5X, т1= 01720. Видно, что размер R0 окрестности источника, в которой задается аналитическое решение, оказывает влияние на величину ошибки: для меньшего R01 ошибка больше, чем для R02, несмотря на меньший шаг по времени.

122

На рисунке 2 показана зависимость от времени специальным образом усредненной по радиусу на кольце 2X < r <3.5X относительной ошибки, которая рассчитывалась по формуле:

p+р

(Sf(tk)) = max X \Fh (rp,p,,tk)-Fa (rp,p,, tk)|

P,

V P=Pi+1

If

Pi+P . . . .

X \Fa (rp,P,,tk)|

V p=Pi+1

Л

(6)

На рис.2а приведена зависимость при t > 200, а на рис.2б в удобном масштабе при t > 3000 показан конец переходного периода и установившийся режим.

В случае проводника (е = е1 = Const, ц = е1 = Const, о = Const > 0) удобно

ввести частоту проводимости Юс=оД SSq ) и толщину скин-слоя

So (®) = у]2/(piPo®a) = а^2/(юю0) .

Дипольное аналитическое решение является точным решением уравнений Максвелла в низкочастотном пределе соа □ со (без учета тока смещения), т. е. является приближенным решением полных уравнений Максвелла (1). Для него магнитное поле определяется формулами (2), а электрическое поле имеет вид:

E (r, P, t ) =

rot H ( r ,P, t)

'(ш-kr)

(k2r 2 -1 - ikr

\ ( X, v) X/

ikr ) v + -—y— (3 + 3ikr - k r

2r 2 )l , (7)

a

причем волновое число является комплексным: & = &(©) = (l-/)/S0(©)eD . Длина волны определяется как X (©) = 2n/Rek (©) = 2тс5о (©). Шаги

пространственной сетки составляли Ar= X/Nx = X/100 и Др = nJNp = п/360 рад = 0.5 град., размер R0 окрестности источника, в которой задается аналитическое решение, был R0 = X. Размеры области моделирования и поглощающего слоя были R1=4.6X, R2=5.4X, т. е. на поглощающий слой также приходится 80 шагов по г. Значение проводимости определялось из условия ©о = 512©. Отметим, что шаг по времени т должен удовлетворять как условию Куранта (4), так и условию аппроксимации затухания т©о ~ 0.1. Были проверены 2 варианта шагов по времени: т1=0/8000 и т2=0.25 т1= 0/32000.

Рис.2. Зависимость от времени рассчитанной по формуле (6) усредненной по радиусу относительной ошибки (8F (t)).

123

Слева на панели а - приведена зависимость при t >200, а справа на панели б - при t >3000 в удобном масштабе показан конец переходного периода и установившийся режим

Рис.3. Тест для проводника при шаге по времени ii=0/8OOO.

Слева а - зависимость от времени рассчитанной по формуле (6) усредненной по радиусу относительной ошибки (bF (t; справа б - зависимость от r рассчитанной по формуле (5) усредненной установившейся относительной ошибки (§F (r, р = 0)) в плоскости поперек диполя (Р=0).

Усреднение проводилось по интервалу [3500; 5000] (в установившемся режиме)

Рис.4. Тест для проводника при шаге по времени т2=0/32000. Обозначения такие же, как и на рис.3

На рис.3 приведены результаты теста для проводника при шаге по времени т1=0/8ООО. На рис.3а показана зависимость от времени рассчитанной

по формуле (6) усредненной по радиусу относительной ошибки (§f (t)), на рис.3б показана зависимость от r рассчитанной по формуле (5) усредненной установившейся относительной ошибки (§f (r,P = 0)) в плоскости поперек диполя (Р=0). Усреднение проводилось по интервалу [3500; 5000]

124

(в установившемся режиме). На рис. 4 приведены аналогичные результаты для в 4 раза меньшего шага по времени т2=0/32ООО. Видно, что уменьшение шага по времени существенно сокращает время установления (и 3000 ^ и 1750), но мало

меняет усредненную установившуюся относительную ошибку (r,P = 0)^.

Заключение

На основании проведенных тестов можно сделать следующие выводы.

1. Окрестность источника, в которой поле задано по формулам аналитического решения, должна быть достаточно большой R0 > 0.57, т. е. содержать ближнюю зону {г < 0.11147} и часть промежуточной (средней) зоны {0.11147, < г < 1.1147}.

2. Дипольное поле в установившемся режиме с хорошей точностью воспроизводится в охватываемой областью моделирования {R0 < г < Rj} части промежуточной (средней) зоны, а также в части дальней зоны {1.1147 < г <3.57}. Нарастание относительной ошибки при {г > 3.57} обусловлено тем, что из-за быстрого убывания по пространству аналитического решения в этой области оно становится очень малым, так что даже при использовании двойной точности в представлениях вещественных чисел даже очень малые ошибки округления сравнимы с точным решением.

Благодарности

Работа выполнена при финансовой поддержке РФФИ, проект № 13-01-00063.

Литература

1. Simpson J. J. Current and future applications of 3-D global Earth-ionospheric models based on the full-vector Maxwell’s equations FDTD method // Surveys Geophys. 2009. Vol. 30. P. 105-130, doi 10.1007/s10712-009-9063-5.

2. Simpson J. J., Taflove A. A review of progress in FDTD Maxwell’s equations modeling of impulsive subionospheric propagation below 300 kHz // IEEE Transactions on Antennas and Propagation. 2007. Vol. 55. No. 6. P. 1582-1590, doi 10.1109/TAP.2007.897138.

3. Две разностные схемы для численного решения уравнений Максвелла для ультра и сверхнизкочастотных сигналов в волноводе Земля - ионосфера

О. И. Ахметов, В. С. Мингалев, И. В. Мингалев и др // Журнал вычислительной математики и математической физики. 2014. Т. 54, № 10. С. 1656-1677.

Сведения об авторах

Ахметов Олег Иршатович,

к.физ.-мат.н., научный сотрудник, Полярный геофизический институт, г. Апатиты, akhmetov@pgia. ru

Мингалев Игорь Викторович,

к.физ.-мат.н., старший научный сотрудник, Полярный геофизический институт, г. Апатиты, [email protected]

Мингалев Олег Викторович,

125

к.физ.-мат.н., и.о. заведующего сектором № 203, Полярный геофизический институт, г. Апатиты, [email protected]

Федоренко Юрий Валентинович,

к.физ.-мат.н., доцент, заведующий сектором, Полярный геофизический институт, г. Апатиты, [email protected]

УДК 537.87, 517.958

М. Н. Мельник, О. И. Ахметов, О. В. Мингалев, И. В. Мингалев

ТОЧНОЕ РЕШЕНИЕ УРАВНЕНИЙ МАКСВЕЛЛА ДЛЯ ТОЧЕЧНОГО ДИПОЛЬНОГО ИСТОЧНИКА В ОДНОРОДНОМ ИЗОТРОПНОМ ПРОВОДНИКЕ И НА ПЛОСКОЙ ГРАНИЦЕ ДВУХ ОДНОРОДНЫХ ИЗОТРОПНЫХ СРЕД

Аннотация

Выведена формула, описывающая решение задачи Коши для телеграфного уравнения в 3-мерном пространстве, аналогичная (и переходящая в нее при нулевой проводимости) формуле Кирхгофа для волнового уравнения. На основе выведенной формулы строится решение задачи о поле электрического диполя Герца с произвольной зависимостью тока от времени в бесконечном однородном изотропном проводнике. Также для случая гармонически зависящих от времени полей получено в наиболее общем случае точное решение задачи о поле горизонтального электрического диполя Герца на плоской границе раздела двух однородных изотропных сред, из которых одна обязательно является проводником.

Ключевые слова:

уравнения Максвелла в проводнике, точные решения, диполь Герца.

M. N. Melnik, O. I. Akhmetov, O. V. Mingalev, I. V. Mingalev

EXACT SOLUTIONS OF MAXWELL EQUATION FOR A POINT DIPOLE SOURCE IN A HOMOGENEOUS ISOTROPIC CONDUCTOR AND ON PLANE BOUNDARY OF TWO HOMOGENEOUS ISOTROPIC MEDIA

Annotation

In this paper a formula has been obtained that describes the solution to the Cauchy problem for the telegraph equation in a three-dimensional space, which is similar to Kirchhoff formula for the wave equation. In the case of zero conductivity this formula transforms to Kirchhoff formula. Based on the derived formula the solution for the problem of the field of the electric Hertz dipole with an arbitrary time dependence of the current in an infinite homogeneous isotropic conductor has been received. Also, for the case of harmonic time-dependent fields an exact solution has been obtained for the problem of the field of the horizontal electric Hertz dipole located on a plane boundary of two homogeneous isotropic media, one of which is necessarily a conductor. Keywords:

Maxwell equations in conductor, exact solutions, Hertz dipole.

Введение

Точные решения уравнений Максвелла для полей от искусственных источников в однородных изотропных средах имеют большое прикладное значение, в том числе для разработки и тестирования методов численного решения уравнений Максвелла. Особенно важны решения

126