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

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

УДК 519.6:551.511

З. В. Суворова, И. В. Мингалев, О. И. Ахметов, О. В. Мингалев

ЯВНАЯ КОНСЕРВАТИВНАЯ СХЕМА РАСЩЕПЛЕНИЯ ДЛЯ УРАВНЕНИЙ МАКСВЕЛЛА

Аннотация

В данной работе представлена новая явная схема для численного интегрирования уравнений Максвелла в изотропных и анизотропных диэлектриках и проводниках. В этой схеме электрическое и магнитное поля вычисляются в одни и те же моменты времени в одинаковых узлах пространственной сетки, а также используется расщепление по пространственным направлениям и физическим процессам. Схема является консервативной, монотонной, имеет 2-й порядок точности по времени и 3-й по пространственным переменным. Представленная схема при моделировании распространения низкочастотных сигналов в волноводе Земля — ионосфера позволяет использовать существенно больший шаг интегрирования по времени, чем широко используемый метод конечных разностей во временной области при одинаковой точности.

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

уравнения Максвелла, схема расщепления, численное моделирование.

Z. V. Suvorova, I. V. Mingalev, O. I Ahmetov, O. V. Mingalev

THE EXPLICIT CONSERVATIVE SPLITTING SCHEME FOR MAXWELL'S EQUATIONS

Abstract

This paper presents a new explicit scheme for the numerical integration of Maxwell's equations in isotropic and anisotropic dielectrics and conductors, as well as in the plasma in the case of the Vlasov-Maxwell system. In this scheme, the electric and magnetic fields are calculated in the same time points in the same spatial grid nodes, and a splitting in spatial directions and physical processes has been used. Scheme is monotonic and has 2nd order accuracy in time and 3rd order accuracy in the spatial variables. The presented scheme allows us to use a much larger step of time integration in modeling the propagation of low-frequency signals in the ionosphere than the widely used method of finite differences in the time domain for the same accuracy. Keywords:

Maxwell's equations, splitting scheme, numerical simulation.

Введение

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

Последние десятилетия для численного интегрирования уравнений Максвелла широко используется метод конечных разностей во временной области (метод КРВО), который в англоязычной литературе принято называть finite-differences time-domain method и использовать сокращение FDTD method [1-6]. Существенным недостатком этого метода является то, что для реализации свободного ухода волн из области моделирования необходимо вводить модельные поглощающие слои толщиной как минимум в десятки шагов сетки, которые прилегают снаружи к границе области моделирования и обеспечивают затухание и торможение сигнала, ушедшего из области моделирования, настолько, чтобы в этой области можно было пренебречь сигналом, отраженным от внешней границы поглощающего слоя. Для определения оптимальных параметров этого слоя, как правило, требуется большой объем тестовых расчетов. Это существенно увеличивает вычислительные затраты, а также усложняет разработку моделей.

На шаг интегрирования по времени в методе КРВО также имеются два ограничения. Первое требует выполнения условия Куранта для обеспечения устойчивости. Второе ограничение связано со схемой учета проводимости среды. Для достижения приемлемой точности требуется, чтобы шаг интегрирования был мал по сравнению с временем изменения полей за счет проводимости среды. При моделировании распространения низкочастотных сигналов в ионосфере (высоты 60-600 км) на пространственной сетке с шагом по высоте 0,5-1 км второе ограничение может быть в 10-100 раз более строгим, чем первое. Это обусловлено тем, что за время прохождения сигналом одного шага сетки поле сигнала существенно затухает и вращается вследствие проводимости ионосферной плазмы и наличия сильного геомагнитного поля, из-за которого ионосферная плазма имеет тензорные диэлектрическую проницаемость и проводимость.

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

Важными достоинствами предлагаемого метода являются консервативность, то есть выполнение в конечно-разностном виде закона сохранения энергии электромагнитного поля и обеспечение равенства нулю дивергенции магнитного поля, а также автоматическое выполнение условий для электромагнитного поля на границе раздела сред, то есть не требуется «сшивки» решения на этих границах.

Построение схемы

Обозначим через г = (л, у, z) декартовы координаты в пространстве К ,

через t - время, через E(r, t) , D(r, t) , H(r, t) и B(r, t) - напряженность и

индукцию электрического и магнитного полей, через j(r, t) - плотность

полного тока в момент времени t в точке r. Будем использовать систему СИ. Запишем уравнение Фарадея и уравнение Максвелла

IB = -rotE(r, t), |D = rot H(r, t) - j(r, t) (1)

в виде гиперболического векторного уравнения 1-го порядка.

Сначала рассмотрим неоднородный изотропный проводник, в котором выполняется закон Ома в форме

j ( r, t ) = ст( r )E( r, t) + J( r, t), (2)

где <r(r) - скалярная проводимость среды, J(r, t) - плотность заданного стороннего тока, а материальные уравнения имеют вид

D(r,t)=£(r)E(r,t), B(r,t)==(r)H(r,t) , (3)

где e(r) и ju( r) — безразмерные относительные диэлектрическая и магнитная проницаемости среды, £ о и ßo — электрическая и магнитная постоянные. Отметим, что в случае диэлектрика или проводника J(r, t) имеет физический смысл плотности стороннего тока в источнике, порождающем сигнал. В случае системы Власова--Максвелла для плазмы а(r) = 0 , £(r) = 1 и ß(r) = 1, а

J(r, t) имеет физический смысл плотности тока плазмы, которая в ходе

численного решения вычисляется по функциям распределения компонент плазмы и на этапе численного решения уравнений Максвелла является заданной функцией.

Обозначим через с0 = l/^/£0ß0 скорость света в вакууме, а через с(r) = с0/yjs(r)/4 r) - скорость света в среде в точке r . Введем перенормированные поля B (r, t) = СоB (r, t) , E (r, t) = £)ß(r) E (r, t) , плотность стороннего тока J ( r, t )= J( r, t )^/l(r) j (soyJs(r ) ) , вектор M(r) = Vc + c(V/(r)^/(r) и функцию fj(r) = ^(r^(££(r)) . Умножая уравнение Фарадея на С0 , получим уравнение:

дВ W Т7Л

— = - rot (сE) .

dt

Подстановка уравнений (2), (3) и введенных обозначений в уравнение Максвелла приводит это уравнение к виду:

d E

— = rot (сВ) -dt

M х В

-^E-J.

(5)

Введем матрицы Rх, R y, Rz и нулевую матрицу O размера 3 X 3 :

Rv

(0 0 0^ 0 0 -1

v0 1 0 у

R

v

( 0 0 1 ^

0 0 0

-1 0 0

R

(0 -1 0 ^ 1 0 0 0 0 0

O

(0 0 0^ 0 0 0

V0 0 0у

Оператор rot в декартовых координатах можно представить в виде

„ da „ da „ Oa

rot a = R x— + R v— + R z—.

dx y dy dz

(6)

(7)

Введем

6-мерные

вектор-столбец

переменных

u = (Bx,By,Bz,Ex,Ey,Ez) и вектор F = (FF2,F3,F4,F5,F6 )T , компоненты

которого заданы формулами:

F = f2 = f3 = 0, (F4,F5,F6f =-

M X В

-VE - J,

(8)

Используя обозначения (6), а также формулы (7) и (8), систему из уравнений (4), (5) можно представить в виде векторного уравнения

du d(AxuVi(avuYd;(azu)

(9)

дt дх^ ' ду ^ у ' д2 где через А х , А у , А 2 обозначены симметричные матрицы размера 6 х 6 , заданные формулами:

А. =

f O cRx^

V cRx O

A y =

( O cRy^ -cRy O

A=

Г O cRz^ cRz O у

(10)

Рассмотрим случай анизотропной среды с тензорными проводимостью а(г) и относительной диэлектрической проницаемостью 8 , но скалярной магнитной проницаемостью:

]{г,() = а{г)Е(г,() + 1(г,(), ), В{г,г) = ^ц{г)Н{г,{).{П)

Поскольку тензор 8 всегда симметричен и имеет положительные собственные значения, то существует любая степень этого тензора, в том числе

-1/ ч - 1/2, ч

тензор 8 \г), обратный к 8 , а также взаимно обратные тензоры 8 (/*] и 8 (г) . В этом случае введем перенормированные поля В (г, / ) = с0В (г, /) ,

плотность

Е (г, I) = 4м(т) г12(Г)Е (г, I) , /(#•,£) =<^//(г) 8/ , тензоры С

стороннего тока

-12

-12

(12)

/\ ^ — 12 / \ — 12 'Щг) = £~ 8 а(г)8 , вектор М(г) = г)//л(г) . С помощью введенных

обозначений и формулы (7) уравнение Фарадея можно записать в виде

д В л

— = -го1; с К с К с К с .

дг ^ > бх^ " / бу^ ' ' " 7

а уравнение Максвелла — в виде

сЬЛ

сК

сК

Г дс с с

(13)

V

-Г| Е-З.

Уравнения (12) и (13) можно записать в виде векторного уравнения (9), в котором симметричные матрицы А х, А у, А 2 определяются формулами:

Л Го сИ

А х=

О КхсЛ

-сИ О

А у=

О Я^с

Ягс

(14)

сИ О

О

а компоненты вектора Р заданы формулами:

ъ = Р2 = ъ = о,

дс с с

-Г\Е-3.

(15)

Векторное уравнение (9) задает 6-мерную линейную гиперболическую систему уравнений 1-го порядка, записанную в консервативной форме. Ее правая часть Р линейно зависит от компонент вектора и . Для численного интегрирования таких систем разработано достаточно много разностных схем, в том числе схемы повышенного порядка точности, которые применяются для уравнений газовой динамики. Наиболее полное описание этих схем содержится в монографиях [7, 8]. Используя метод расщепления по пространственным направлениям и физическим процессам [8-10], можно построить явные монотонные схемы численного интегрирования системы (9), которые сводятся к последовательному интегрированию 1 -мерных по пространству гиперболических систем уравнений.

Для изотропной среды и для анизотропной среды с симметричным тензором проводимости (при отсутствии холловской проводимости) на каждом временном шаге нужно последовательно проинтегрировать 3 системы уравнений, например в таком порядке:

^ + —(Ахи') = Рх , ^ + ^(Ауи») = Ру , ^ + — (А,и''') = Р, . (16)

Ы дх\ х > х Ы ду\ у / у Ы дЛ 2 ) 2

При этом правые части систем (16) выбираются так, чтобы они не изменялись на своем шаге расщепления и удовлетворяли равенству

\

Ж = Жх + Жу + . (17)

Оптимальный выбор этих правых частей будет обсуждаться ниже. На каждом из шагов расщепления для двух компонент магнитного поля и двух компонент электрического поля, ортогональных направлению шага, рассчитывается распространение сигнала конечно-разностным способом, а также рассчитывается по аналитическим формулам затухание за счет проводимости третьей компоненты электрического поля. В качестве начальных условий для каждой системы уравнений в (16) берутся значения, рассчитанные в результате предыдущего шага расщепления. Сохранить второй порядок аппроксимации по времени в схеме расщепления можно, если циклически изменять порядок выполнения шагов расщепления. Например, выполняя сначала в следующем порядке шаги по пространственным направлениям: Ху2 , уХ2 , 2Ху , Х2у , у2Х , 2уХ . Обоснование этого утверждения содержится, например, в монографиях [8, 9].

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

с)Е т

—=[лхе], а={1у2,,1ху) ,

в которой , , 71у2 являются компонентами антисимметричной части

тензора ^ . Эта система задает вращение поля Е с вектором угловой скорости

&. Соответствующее циклическое изменение последовательности выполнения шагов расщепления обеспечивает второй порядок аппроксимации по времени и в этом случае.

Пусть заданы равномерная сетка по времени tn = ¿0 + пт, где Т - шаг по

времени, и равномерная пространственная сетка в декартовых координатах, целые и полуцелые узлы которых заданы соотношениями

Х = Х0 + 1 ■ НХ. у.1 = уо + ] ■ Ну , 2к = 20 + к ■ Гг],к = (Х. ^^. 2к ) ,

Х1+1/2 = Х0 + (1 + 1 / 2) ■ К . у]1/2 = уо + О + 1 / 2) ■ Ну .

*к+1/2 = 20 + (к + 1/2) ■ К,

в которых КХ, Ну , Н2 - шаги сетки по осям Х, у, 2 соответственно. Будем использовать для значений функции / в узлах сетки обозначения

/ (Гг,],к, tn ) = /и],к .

Рассмотрим для случая изотропной среды послойный переход для первой системы из (16), которая интегрируется на шаге расщепления по направлению х . В этом случае представим заданную первой формулой в (10)

-1

матрицу А х в виде Ах = ОЛQ • с , где Л — диагональная матрица, на диагонали которой стоят поделенные на с собственные числа матрицы

О — матрица, столбцы которой есть правые собственные векторы матрицы

-1

, определенные с точностью до множителя, а О есть матрица, обратная

-1

к О . Строки матрицы О есть левые собственные векторы матрицы А х . Эти матрицы для изотропной среды постоянны. Их можно взять в виде

Л =

(0 0 0 0 0 0^

0 1 0 0 0 0

0 0 -1 0 0 0

0 0 0 0 0 0

0 0 0 0 1 0

ч0 0 0 0 0 -1,

Q =

(1 0 0 0 0 0 ^

0 1/2 0 0 0 1/2

0 0 1/2 0 1/2 0

0 0 0 1 0 0

0 0 -1/2 0 1/2 0

ч0 -1/2 0 0 0 1/2,

О =

(1 0 0 0 0 0^

0 1 0 0 0 -1

0 0 1 0 -1 0

0 0 0 1 0 0

0 0 1 0 1 0

у0 1 0 0 0 1,

С помощью введенных матриц Л , О и О первую систему уравнений в (16) представим в виде векторного уравнения ды д ( -1 \

Введем вектор-столбец характеристических переменных (линеаризованные инварианты Римана) по направлению х , заданный формулой

-1

№ = О ы . Компоненты вектора № заданы формулами:

= Вх , ^2 =Ву - Ег , = Вг - Е

г у ■

= Ех , = Вг + Еу , = Ву + Ег

(19)

-1

Умножим уравнение (18) слева на матрицу О . В результате получим векторное уравнение

дж А д( сж) 1

—+л——-=о

д дх

(20)

которое является системой из шести скалярных уравнений для характеристических переменных. Отметим, что для изотропной среды правые части систем (16) лучше всего взять в виде

^ = (0,0,0,-ЛЕХ - Зх, -Ы2Бх, ЫуБх) Ру = (0, о, 0, МгБу,-ЛЕу - Зу, -МхБу )Т,

Т

Р = (0,0,0,-МуБг, МХБ2,-ЛЕг -)

что обеспечивает выполнение формул (17) и (8), а также наиболее простую форму уравнений для характеристических переменных на каждом шаге расщепления. Для изотропной среды уравнения системы (20) с учетом формул (17) и (19) имеют вид:

-71 = 0, —4 = -ф*-4 - Зх , (22)

дt дt

= -муЩ, + дА™А = -мл, (23)

дt дх дt дх

—^--^—^ = М271, —^--^—^ = Мы , (24)

дt дх дt дх

Отметим, что в силу первого уравнения в (22) на шаге расщепления

по X постоянна 7 1 и, следовательно, постоянны правые части уравнений (23),

(24). Также отметим, что второе уравнение в (22) является линейным уравнением 1-го порядка с постоянным коэффициентом и его решение определяется по формуле Коши:

т

<+1= < ехр(-^т)-^ Jx(tn + я) ехр(-ц(т- я)) ds , (25)

0

в которой интеграл можно вычислять численно или аналитически, если задана аналитическая зависимость стороннего тока J(г, t) от времени. Эта формула

определяет послойный переход для переменной 74 . В итоге для изотропной

среды численное интегрирование первой системы из (16) сводится к численному интегрированию четырех независимых 1 -мерных по пространству уравнений переноса (23) и (24).

Уравнения (23) описывают волны, бегущие вдоль оси X слева направо, а уравнения (24) описывают волны, бегущие вдоль оси X справа налево. Для этих уравнений предложено большое число разностных схем [7, 8]. Существует несколько явных монотонных схем, которые имеют 3-й порядок точности по пространству и 1 -й порядок точности по времени. Тестовые расчеты показали, что для получения приемлемого качества численного решения необходимо использовать схемы повышенного порядка точности. Авторами данной работы была предложена следующая явная гибридная схема, которая имеет 2-й порядок точности по времени и 3-й по пространству.

Рассмотрим построение этой схемы для уравнений (23). Эти уравнения можно записать в виде

а/ а( е/)

а +лаег=£• (26)

где правая часть £ не зависит от времени. Явные разностные схемы для уравнения (26) записываются в виде

/1]к - /1]к = - (С+1/2]к / +1/2]к - С1 -1/2]к /1-1/2]к ) + Т£1]к • (27) где С1+1/2]к = е1+1/2]кт1 К > 0 есть число Куранта в узле сетки Г+1/2,},к ■

Порядок точности схемы (27) и ее монотонность обеспечиваются способом вычисления значений функции / в полуцелых по X узлах сетки. Для записи этого способа обычно используют параметр, называемый анализатором гладкости функции / и заданный формулой

гП гП

Ж = 1цк - Л-ик •

гП _ гП

JI +1 ]к Jук

и функцию 0(Ж) , называемую ограничителем потока. Формулы для

гП

вычисления /¿+1/2 записывают в виде

-я+2(>-с+*к) ^ (/+„к-/ик) (29)

Для схемы Лакса — Вендроффа О(Ж) = 1, для схемы Бима — Уорминга Ж) = Ж . Эти схемы имеют 2-й порядок точности по времени и по пространству. При О(Ж) = 0 получается монотонная схема 1-го порядка точности. Как показано в работах [11, 12], для монотонности схемы (27), (29) необходимо выполнение следующего условия на ограничитель потока:

0 < О(Ж) < Ш1п(2, Ж+1Ж |). (30)

Из (30) следует, что схема Лакса Вендроффа монотонна при выполнении условия Ж>0.5 , а схема Бима — Уорминга монотонна при выполнении условия 0 < Ж < 2 . Эти схемы являются немонотонными, однако области монотонности этих схем перекрываются. В работе [13] была предложена явная гибридная схема для численного интегрирования уравнений динамики несжимаемой жидкости, в которой осуществляется переключение между схемами Лакса — Вендроффа и Бима — Уорминга. В предложенной авторами гибридной схеме осуществляется переключение между шестью схемами, которые имеют 2-й порядок точности по времени и 3-й по пространству и которые монотонны при выполнении условий переключения. Эти условия, а также сама схема записываются в следующем виде:

1. Если f"+lJk — fyk = 0, то полагаем О(ЭТ) = 0 .

2. Если f»+ljk — j ф 0 и ЭТ< 0, то полагаем в(ЭТ) = (l — C1+m jk )/l.

3. Если tfj — j Ф0 и ЭТ>0, то:

а) при 0 <ЭТ<(1 — Ci+mjk )/(3 — Ci+mjk) полагаем О(ЭТ) = ЭТ(2 — CMI2jk ),

б) при (l — C+1/2 ]к )/(3 — C+i/2 ]к ) < ЭТ < (2 — C+1/2 ]к )/(5 — C+i/2 ]к ) полагаем О(ЭТ) = ЭТ (l + Ci+1/2]к)/2 + (l — C+1/2д )/2,

в) при (2 — Ci +l/2 jk )/ (5 — Ci+l/2 jk ) < ЭТ < (4 + Ci+l/2 jk )/ (1 + Ci+l/2 jk ) полагаем О(ЭТ) = ЭТ(1 + Ci+1/2 k)/3 + (2 Ci+1/^.,-k V3,

г) при (4 + CMI2jk)/(l + CMI2jk)<ЭТ полагаем С(ЭТ) = (3 + C1+mß)/2.

Отметим, что пункты а) и б) являются уточнениями схемы Бима — Уорминга, а пункт г) является уточнением схемы Лакса — Вендроффа до 3-го порядка точности по пространству. Пункт в) является схемой 3-го порядка точности по пространству, использующей центральные разности. Прямой проверкой можно убедиться, что условие монотонности (30) выполняется для пунктов а), б), в) и г).

Рассмотрим уравнения (24) и запишем их в виде

df d( cf)

f — g' (31)

dt dx

где правая часть gi не зависит от времени. Явные разностные схемы для уравнения (31) записываются в виде

fijk+ — fijk = (Ci+1/2jk fi+l/2jk — Ci —1/2jk fi—l/2jk ) + Tgijk > (32)

который аналогичен (27). В схеме (32) значения функции f в полуцелых по X узлах сетки вычисляются по формуле

fi+l/2 jk = fi+ljk — 1 (1 — Ci+1/2 jk ) &(ЭТ) (fi+l jk — fFjk ) '

в которой анализатор гладкости функции f задан формулой

гП _ гП

ЭТ = ——j-i—— , а условия переключения для ограничителя потока G(ЭТ)

гП _ гП

Ji+ljk Jijk

точно такие же, как в изложенной выше схеме для уравнения (26).

Рассмотрим задание граничных условий для предложенной схемы на шаге расщепления по направлению X . Поскольку уравнения (23) описывают волны, бегущие вдоль оси x слева направо, то необходимо задавать значения характеристических переменных w и W5 на левой границе. Аналогично,

поскольку уравнения (24) описывают волны, бегущие вдоль оси x справа налево, то необходимо задавать значения характеристических переменных w3 и w6 на правой границе.

Условие свободного ухода волн на левой границе области моделирования реализуется в случае, когда на этой границе задаются значения

переменной W2 , не связанные со значениями переменной w6 , и значения переменной W3 , не связанные со значениями переменной w5 . Условие свободного ухода волн на правой границе области моделирования реализуется в случае, когда на этой границе задаются значения переменной w5, не связанные

со значениями переменной W3 , и значения переменной w6 , не связанные со значениями переменной W2 . Условие отражения волн на левой границе области моделирования реализуется, когда на этой границе задаются значения переменных W2 = KjWe и W3 = KjW5 , где K - коэффициент отражения на

левой границе. Аналогично условие отражения волн на правой границе области моделирования реализуется, когда на этой границе задаются значения переменных w6 = Krw2 и w5 = KrW3 , где Kr - коэффициент отражения на правой границе.

Послойные переходы для второй системы из (16), которая интегрируется на шаге расщепления по направлению y , и для третьей системы из (16), которая интегрируется на шаге расщепления по направлению z , выполняются аналогично описанному переходу на шаге расщепления по направлению x .

В случае анизотропной среды на шаге расщепления по направлению x -1

матрицы Q и Q имеют несколько более сложный вид, чем в случае изотропной среды, но уравнения системы (18) полностью аналогичны таковым в случае изотропной среды. Можно показать, что в случае анизотропной среды

W = Bx , выполняется первое из уравнений (20), w4 является линейной комбинацией компонент электрического поля, уравнение для w не содержит

производную по x , а остальные четыре уравнения являются уравнениями переноса, аналогичными уравнениям (21)-(24).

Заключение

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

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

Проведенные авторами тестовые расчеты показали монотонность и консервативность предложенной схемы. Для нее представлены граничные условия свободного ухода волны из области моделирования. Представленная схема при моделировании распространения низкочастотных сигналов в волноводе Земля — ионосфера позволяет использовать существенно больший шаг интегрирования по времени, чем широко используемый метод конечных разностей во временной области при одинаковой точности расчетов.

Благодарности. Работа выполнена при финансовой поддержке РФФИ, проект 17-01-00100.

Литература

1. Numerical solution of initial boundary value problems involving Maxwell's equations in isotropic media / Yee Kane // IEEE Transactions on Antennas and Propagation. 1966. Vol. 14. P. 302-307.

2. 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.

3. 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 (June 2007). P. 1582-1590.

4. Paul D. L., Railton C. J. Spherical ADI FDTD method with application to propagation in the Earth ionosphere cavity // IEEE Transactions on Antennas and Propagation. 2012. Vol. 60, No. 1 (January 2012). P. 310-317.

5. Yu Y., Simpson J. J. An collocated 3-D FDTD model of electromagnetic wave propogation in magnetized cold plasma // IEEE Transactions on Antennas and Propagation. 2010. Vol. 58, No. 2 (February 2010). P. 469-478.

6. Семенов А. Н., Смирнов А. П. Численное моделирование уравнений Максвелла с дисперсными материалами // Математическое моделирование. 2013. Т. 25, № 12. С. 19-32.

7. Semenov A. N., Smirnov A. P. Numerical modeling of Maxwells equations with dispersive materials // Matem. Mod. 2013. Vol. 25, No. 12. S. 19-32.

8. Куликовский А. Г., Погорелов Н. В., Семенов А. Ю. Математические вопросы численного решения гиперболических систем уравнений. 2-е изд., испр. и доп. М.: Физматлит, 2012. 656 с.

9. Бисикало Д. В., Жилкин А. Г., Боярчук А. А. Газодинамика тесных двойных звезд. М.: Физматлит, 2013. 632 с.

10. Ковеня В. М., Яненко Н. Н. Метод расщепления в задачах газовой динамики. Новосибирск: Наука, 1981.

11. Вязников К. В., Тишкин В. Ф., Фаворский А. П. Построение монотонных разностных схем повышенного порядка аппроксимации для систем уравнений гиперболического типа // Матем. Моделирование. 1989. Т. 1, № 5. С. 95-120.

12. Harten A. High resolution schemes for hyperbolic conservation laws // J. Comp. Phys. 1983. Vol. 49. P. 357.

13. Sweby P. K. High Resolution Schemes Using Flux Limiters for Hyperbolic Conservation Laws // SIAM J. Numer. Anal. 1984. Vol. 21. P. 995.

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

Суворова Зоя Викторовна

программист, Полярный геофизический институт, Апатиты

E-mail: suvorova@pgia.ru

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

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

E-mail: mingalev_i@pgia.ru

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

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

E-mail: mingalev_o@pgia.ru

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

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

E-mail: akhmetov@pgia.ru