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

УДК 519.6: 551.511

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

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

Аннотация

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

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

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

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

TESTING THE EXPLICIT CONSERVATIVE SPLITTING SCHEME FOR MAXWELL'S EQUATIONS

Abstract

This paper presents the test results of the new explicit scheme for numerical integration of Maxwell's equations in isotropic and anisotropic dielectrics and conductors. With the help of test calculations proved the monotonicity and conservatism of the scheme, the dependence of the angle of incidence of the wave reflection coefficient of the wave from the border region modeling. It was also made the comparison of the numerical solution with analytical solution of the problem about a field of a point electric dipole on a flat boundary between two media, as well as problem on the field of an electric dipole in a homogeneous conductor. Keywords:

Maxwell's equations, splitting scheme, numerical integration. Введение

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

Предложенная схема и условия свободного ухода волны на границе области моделирования без использования модельного поглощающего слоя дали следующую зависимость от угла падения отношения амплитуды волны, отраженной от границы области моделирования, к амплитуде падающей волны. Для волн, падающих под углом от 80° до 90°, указанное отношение не превышает 0,01, при угле падения 60° это отношение уже составляет примерно 0,05, при угле падения 45° — примерно 0,16, при угле падения 27° — примерно 0,35, а при угле падения 18,4° — примерно 0,43.

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

Монотонность схемы

Монотонность схемы была проверена авторами при моделировании распространения плоских и сферических волн с П-образным профилем. Область моделирования представляла собой куб, шаги равномерной пространственной сетки были одинаковы по всем направлениям. Число узлов сетки по каждому направлению было равным 385. На границе области моделирования было задано условие свободного ухода волны. Для представления численного и аналитического решений используются безразмерные координаты и безразмерные компоненты напряженности электрического поля. На рис. 1 приведены зависимости от координаты X компоненты поля Еу в бегущей вдоль оси X слева направо

плоской линейно поляризованной волне с П-образным профилем в момент времени, когда передний фронт волны еще не достиг границы области моделирования. Шаг пространственной сетки по каждому направлению был равным 0,025, а число Куранта было равным 0,25. Пунктирная линия представляет численное решение, а сплошная линия — аналитическое решение. Из рис. 1 видно, что численное решение хорошо воспроизводит амплитуду, передний и задний фронты волны, что отсутствуют нефизические осцилляции численного решения и что скачок в численном решении размазывается не более чем на 8 узлов сетки. Расчеты показали, что эта ширина размазывания скачка сохраняется во времени.

Рис. 1. Безразмерная компонента поля Еу в плоской линейно поляризованной волне с П-образным профилем, бегущей вдоль оси X слева направо: 1 — точное решение; 2 — численное решение

Консервативность схемы

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

Также в этих тестовых расчетах проверялся баланс энергии для контрольного объема каждого узла сетки и баланс энергии для всей области моделирования. В безразмерном виде в целых узлах сетки в целые моменты

времени вычислялись энергия поля и = Е2+ Б2 и вектор потока энергии

Л = [Е х Б] , а также Л1уЛ с помощью центральных разностей. Уравнение

локального баланса энергии при отсутствии источников излучения записывается

ди .. Л

в виде -+ шуЛ — 0 . Относительное нарушение этого уравнения в целых узлах

д/

сетки вычислялось в полуцелые моменты времени по формуле

&

п+1/2 _

2 (ип+1 - ип) д+л у (лп + Лп+1 У (2 \ип+1 - ип\ /т + |Лу (лп + Л"+1 )|)

для тех узлов сетки, в которых 2|ип+1 -ип|Т + |&у(Лп + Лп+1 0 . При

вычислениях с одинарной точностью величина &п+12 не превышала 10-5 во все моменты времени.

Для проверки общего баланса энергии в безразмерном виде в целые моменты времени с двойной точностью вычислялись суммарная энергия поля во всей области моделирования Ш и суммарный поток энергии поля через границу области моделирования Р . При этом поля, энергия поля и ее поток в узлах сетки вычислялись с одинарной точностью. Уравнение общего баланса энергии в области моделирования при отсутствии источников излучения записывается

дШ п_п Л

в виде -+ Р — 0 . Относительное нарушение этого уравнения вычислялось

д/

в полуцелые моменты времени по формулам

8п+12 = |2 (Шп+1 - Шп )/т+ Рп + Рп+1|/ (2\жп+1 - Т + \Рп + Рп+11).

Расчеты показали, что величина Зп+12 не превышала 10-6 во все моменты времени. Таким образом, тестовые расчеты показали, что предложенная схема обеспечивает выполнение закона сохранения энергии с высокой точностью.

Численное отражение волн от границы области моделирования

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

для волн, падающих под углом 80-90°, отношение не превышает 0,01;

при угле падения 60° это отношение уже составляет ~ 0,05;

при угле падения 45° — примерно 0,16;

при угле падения 27° — примерно 0,33;

при угле падения 18,4° — примерно 0,43.

На рис. 2 в безразмерном виде представлены результаты расчетов падения плоской линейно поляризованной волны на границу области моделирования под углами 63°, 45° и 27°. Электрическое поле волны направлено вдоль оси X . Падающая волна движется в плоскости у2 в сторону

левого нижнего угла области моделирования. Безразмерная амплитуда изменения электрического поля в этой волне равна 1. На левой части рисунка волна падает под углом 45°, а на правой части она падает под углом 63° к левой по у границе и под углом 27° к нижней по 2 границе. Видно, что на левой части рис. 2 амплитуда отраженной волны не превышает 16 % от амплитуды падающей волны. Также на правой части рис. 2 видно, что для угла падения 27° амплитуда отраженной волны не превышает 35 % от амплитуды падающей волны, а для угла падения 63° — 6 %.

__У

^_:_I_: •:"":" к 1

-0.2 0 0.2 0.4 0.6 0.8 1 -0.3 -0.1 0 0.2 0.4 0.6 0.8 1 Рис. 2. Безразмерное электрическое поле в плоской волне, линейно поляризованной вдоль оси x и падающей на границу области моделирования в плоскости у2 , и в волне, отраженной от этой границы. Справа угол падения равен 45°, слева этот угол равен 63° на левой по у -границе и равен 27°

на нижней по 2 -границе

Расчеты показали, что для полного подавления эффекта отражения волн, падающих почти по касательной к границе, достаточно ввести модельный поглощающий слой толщиной 25-30 узлов.

Сравнение результатов расчетов с точными решениями

Авторами были получены численные решения задач о распространении сигнала от точечного электрического диполя, расположенного на границе однородных диэлектрика и проводника, а также расположенного в однородных проводнике и диэлектрике. Эти численные решения сравнивались с точными аналитическими решениями этих задач.

Область моделирования представляла собой куб, шаги равномерной пространственной сетки были одинаковы по всем направлениям, диполь и начало координат находились в центре области моделирования. Направление электрического диполя было постоянным и совпадало с направлением оси X . Дипольный момент менялся во времени по формуле гармонических колебаний с амплитудой, которая в начальный момент времени равнялась нулю, а затем в течение двух периодов колебаний плавно возрастала от нуля до некоторого значения и далее оставалась постоянной. Таким образом, дипольный излучатель плавно включался в начальный момент времени и плавно выходил на режим гармонических колебаний с постоянной амплитудой. Поперечный размер области моделирования составлял 8 длин волн установившихся колебаний в диэлектрике. На границе области моделирования было задано условие свободного ухода волны. Число узлов сетки по каждому направлению было равным 385 (49 узлов на длину волны), а число Куранта было равным 0,5. Для представления численного и аналитического решений используются безразмерные координаты X, у, z , выраженные в длине волны установившихся колебаний, безразмерные компоненты индукции магнитного поля. В случае двух сред верхняя половина области моделирования (при z > 0) являлась диэлектриком, а нижняя половина (при z < 0) являлась проводником.

Задача о поле точечного электрического диполя, расположенного на границе двух однородных сред, в случае, когда дипольный момент гармонически

колеблется с амплитудой Рт и частотой т0 , имеет точное аналитическое

решение [1, 2]. В этом решении векторный потенциал магнитного поля А(г, ?)

имеет только две ненулевые компоненты вдоль осей x и z . Их комплексные амплитуды в системе СИ заданы формулами

и и и* Р +ГЛе~' ^^ J0 Х^ТУ?) йХ

иош*0Рт Г--у_)_ А (г) = 0,

0 _ (1)

= и иит РтХ z 1кии4X2 + у2 Л Е(Х)(к2иЛХ2 - к2 + к22иАХ2 - к2 ) '

в которых Е(Х) = и ^Х — к22 + и2 VХХ — к,2 , ^ (^) и ^ (^) — функции

Бесселя, индекс ( = 1 используется для верхнего полупространства ^ > 0} , в котором магнитная и диэлектрическая проницаемости среды и ее проводимость равны соответственно , квадрат волнового числа

к1 = ^^„^о ^^о^о — /г) , а индекс а = 2 используется для нижнего полупространства {г < 0} , в котором указанные параметры среды равны ¿2,^2'Г2 , = — Т)• Подставив формулы (1) в формулу

В (г, г) =го1 А( г, г) , получим квадратурную формулу для комплексной амплитуды компоненты поля В . Из-за больших размеров эту формулу не

приводим. Расчеты по этой формуле проводились с помощью численного вычисления входящих в нее интегралов.

На рис. 3 представлена зависимость безразмерной амплитуды установившихся колебаний компоненты поля В от безразмерной координаты X для численного решения задачи о поле гармонически колеблющегося точечного электрического диполя на границе диэлектрика с параметрами /Л1 =е1 = 1, Г1 = 0 и проводника с параметрами ¿2=е2 = 1, г2 = е0ю0, а также эта зависимость,

рассчитанная по квадратурной формуле для точного аналитического решения вдоль двух прямых, проходящих параллельно оси X через точки с безразмерными координатами у = 0,5, г = —1 в проводнике и у = 0,5, г = 1 в диэлектрике. Видно, что две указанные зависимости очень близки между собой. Этот факт подтверждает хорошую точность предложенной схемы.

Рис. 3. Амплитуда гармонических колебаний компоненты поля В

при наличии точечного диполя на границе двух сред. Слева в проводнике, справа в диэлектрике, 1 — точное решение, 2 — численное решение

В однородном проводнике векторный потенциал точечного электрического диполя, находящегося в начале координат и направленного вдоль оси х , в системе СИ задан формулой

( х,' ) = —

¿¿0 4^1 х

(

ехр

— х| 1 йР (г — \х\/с )

2 с

йг

+

+

а — | х[

Г I ( — (Iх|/с)2 1 еХР (~\Т'2) йР (г — Т) йТ

в которой использованы обозначения: 2 = ст/(££0) , где ст — проводимость среды, с = \Це0ЛоЛ — скорость света в среде, через Р(г) обозначен момент диполя как функция времени, а через /0 (5) =М0 (¿5), ^ (5) = —¿М1 (¿5)

— модифицированные функции Бесселя. В случае, когда диполь направлен вдоль оси х, для компонент индукции магнитного поля вытекают следующие формулы:

Бх(х, г) = 0, ву(X, г) = Г Ф(|х\, г), Бг (х, г) = —у Ф(|х|, г),

Ф(| х|, г)

ЛЛо2

Ъ2же 2 |х

-ехр I

щ х

2с

ёР (г — |х|/с)

ёг

4жс\ х'

> 2с Л аг2 4 У V

л

лл0 ( щ\х\ V а2 р (г— И/с) Г с ар (г—|х|/с)

2 ехр--^ --— ' + — + — —-— 7

2с /л |х| 2 ...

Ч I У У

2

г

ёР1г—£) ехр Г—21'

\б^с J ёг V 2

ЛЛо2

I

х1/с

,г

^Щ^т2 —(| х|/ с )2

(2)

ёт.

-4 -3.5 -3 -2.5 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 2.5 3 3.5 4

z, длина волны

Рис. 4. Зависимость выражения (х2 + у2 + г2)Бу от г при х = 0 и у = 0,5

в момент времени, когда передний фронт сигнала находится в области

моделирования

На рис. 4 представлены численное решение задачи о поле точечного электрического диполя в однородном проводнике и точное аналитическое решение, рассчитанное по квадратурным формулам с помощью численного

вычисления интегралов. Приведена зависимость произведения (X2 + у2 + г2)В,

в котором все величины безразмерные, от безразмерной координаты 2 при значениях безразмерных координат X = 0 и у = 0,5 в момент времени, когда передний фронт сигнала еще не достиг границы области моделирования. Пунктирная линия представляет численное решение, сплошная линия - точное аналитическое решение, рассчитанное по формулам (2), а штрихпунктирная линия - разность этих решений. Видно, что численное решение хорошо воспроизводит передний фронт сигнала и что точки максимумов и минимумов

произведения (X2 + у2 + г2) В для численного и аналитического решений

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

(х2 + у2 + г2^В по переменной 2 имеет максимальные значения.

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

Заключение

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

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

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

Литература

1. Тихонов А. Н., Самарский А. А. Уравнения математической физики. 5-е изд-е. М.: Наука, 1977.

2. Светов Б. С. Основы геоэлектрики. М.: Изд-во ЛКИ, 2008.

Сведения об авторах Суворова Зоя Викторовна

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

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

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

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

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

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

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