CC BY
23
УДК 537.87, 517.958
З. В. Суворова, И. В. Мингалёв, О. В. Мингалёв, О. И. Ахметов ЯВНАЯ СХЕМА РАСЩЕПЛЕНИЯ ДЛЯ УРАВНЕНИЙ МАКСВЕЛЛА Аннотация
В данной работе представлено новое семейство явных схем для численного интегрирования уравнений Максвелла в диэлектриках, проводниках с законом Ома, а также в сильно- и слабоионизированном газе. В этих схемах электрическое и магнитное поля вычисляются в одни и те же моменты времени в одинаковых узлах пространственной сетки, а также используется расщепление по пространственным направлениям и физическим процессам. Схемы имеют 2-й порядок точности по времени и 3-й — по пространственным переменным и являются монотонными. Также они позволяют реализовать большой набор граничных условий, в том числе свободный уход волны без использования модельных поглощающих слоев.
Ключевые слова:
уравнения Максвелла, численное интегрирование, схема расщепления.
Z. V. Suvorova, I. V. Mingalev, O. V. Mingalev, O. I. Akhmetov THE EXPLICIT SPLITTING SCHEME FOR MAXWELL'S EQUATIONS Abstract
This paper presents a new family of explicit schemes for the numerical integration of Maxwell's equations in dielectrics and in conductors, Ohm's law, as well as in weakly and highly ionized gas. In these schemes the electric and magnetic fields are calculated in the same time points in the same spatial grid points, and uses a splitting in spatial directions and physical processes. Schemes have 2nd order accuracy in time and 3rd order accuracy in the spatial variables and are monotonic. They also allow you to implement a large set of boundary conditions, including free care wave model without the use of absorbing layers.
Keywords:
Maxwell's equations, numerical integration, splitting scheme. Введение
Существует широкий круг задач, в которых требуется проводить численное интегрирование уравнений Максвелла. К таким задачам относятся: численное моделирование динамики бесстолкновительной плазмы в рамках системы уравнений Власова — Максвелла, распространение сигналов в диэлектриках и проводниках, в которых выполняется закон Ома, а также в разреженном слабоионизированном газе, например в ионосфере Земли.
Последние десятилетия для численного интегрирования уравнений Максвелла широко используется метод конечных разностей во временной области (метод КРВО), который в англоязычной литературе принято называть finite-differences time-domain method и использовать сокращение FDTD method. Описание этого метода приведено в работах [1-6]. Он использует представление
уравнений Максвелла в интегральной форме по пространственным переменным и дифференциальной по времени, а также разностную сетку из ячеек Йее (Yee lattice) (см. [1-6]). Этот метод относительно прост в реализации и позволяет строить для уравнений Максвелла явные двухслойные по времени разностные схемы, которые имеют 2-й порядок точности как по времени, так и по пространству (см., например, обзоры [2, 3] и статьи [4-6]) и устойчивы при выполнении условия Куранта. Важным достоинством метода является автоматическое выполнение условий для электромагнитного поля на границе раздела сред, то есть не требуется сшивки на этих границах.
В данном методе используется регулярная пространственная сетка в декартовых координатах. Его главные особенности заключаются в том, что магнитное поле вычисляется в моменты времени, смещенные на полшага интегрирования относительно моментов времени, в которые вычисляется электрическое поле, и в том, что узлы пространственной сетки, в которых вычисляется магнитное поле, смещены на полшага относительно узлов, в которых вычисляется электрическое поле. Для аппроксимации пространственных производных используются центральные разности.
Существенный недостаток метода заключается в сложности реализации граничного условия свободного ухода волн из области моделирования, что требуется во многих прикладных задачах. Для реализации такого граничного условия вводят модельные поглощающие слои, которые обеспечивают затухание отраженного от границы области моделирования сигнала. Наличие этих слоев значительно увеличивает вычислительные затраты и усложняет разработку моделей.
В данной работе предложен новый метод численного интегрирования уравнений Максвелла, в котором электрическое и магнитное поля вычисляются в одни и те же моменты времени в одних и тех же узлах пространственной сетки, а также используется противопотоковая аппроксимация пространственных производных и схема расщепления по пространственным направлениям и по физическим процессам. Этот метод имеет второй порядок точности по времени, третий по пространству и допускает реализацию любых граничных условий без модельных поглощающих слоев.
Построение схемы в проводнике
Рассмотрим уравнения Максвелла в системе СИ в декартовых координатах г (х, у, ¿) в неоднородном проводнике, в котором выполняется закон Ома. Обозначим через Е(г) напряженность электрического поля, через В(г) — индукцию магнитного поля, через е(г) и ц(г) — диэлектрическую и магнитную проницаемости среды, через е0 и ц — эти же проницаемости вакуума, через
с(г) = . — — скорость света в проводнике в точке г, а через о(г) —
фо^Лг )р(г )
скалярную проводимость среды. Введем перенормированную индукцию
магнитного поля B(r) = c(r)B(r) и вектор M(r) = -
c(r) (Vs V^
2
s ^
. С помощью
введенных обозначений уравнение Фарадея и уравнение Максвелла можно записать в следующем виде:
дБ „ дБ
— = -с гол Б, —
дл д
= с гол Б +
М х Б
Б .
(1)
Введем 6-мерные вектор-столбец переменных и = ( Вх, Ву, Вг, Ех, Еу, Ег )
векторы потоков вдоль координатных осей
X (и) = (>
= 10, -Е., Е. ,0, В2, -В
Г (и) = (Е2,0, -Ех, -В, ,0, Вх ) , Z (и) = (-Еу, Ех ,0, Ву, -Вх ,о)Т , а также вектор
правых частей Б (и) = ( К, К, К, К, К, К ) , компоненты которого заданы формулами
К = К = К = 0.
К = МуВ, -Ы2ВУ -
стЕ„
К = МВХ -МВ2 -■
(гЕ„
К = МВу -МВХ -
аЕ„
5 = М2ВХ МХВ2 , = МхВу
Используя эти векторы, уравнения (1) можно записать в виде векторного уравнения:
ди дХ дГ дZ ^
— + с— + с— + с— = Б . дл дх ду дг
(2)
Отметим, что векторы потоков можно представить в виде X(и) = Ахи,
Г (и) = Ауи, Z(u) = А2и, где Ах =
дХ (и )
Ау =
дГ(и) , дZ (и)
А2 = —есть
ди ди ди
постоянные симметричные матрицы Якоби векторов потоков
(0 0 0 0 0 0 ^ (0 0 0 0 0 1 ^ (0 0 0 0 -1 0^
Ах =
0 0 0 0 0 -1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 -1 0 0 0 0
Ау =
у
0 0 0 0 0 0 0 0 0 -1 0 0 0 0 -1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0
Аг =
у
0 0 0 1 0 0
0 0 0 0 0 0
0 1 0 0 0 0 -1 0 0 0 0 0
0 0 0 0 0 0
у
Векторное уравнение (1) задает 6-мерную гиперболическую систему уравнений 1-го порядка. Для численного интегрирования таких систем разработано достаточно много разностных схем, в том числе схемы повышенного порядка точности, которые применяются для уравнений газовой динамики. Наиболее полное описание этих схем содержится в монографиях [7, 8]. Используя метод расщепления по пространственным направлениям [8, 9], можно построить явные монотонные схемы численного решения системы (1), в которых численное интегрирование этой системы сводится к последовательному интегрированию одномерных по пространству гиперболических систем
уравнений. На каждом временном шаге нужно последовательно проинтегрировать 3 системы уравнений. Вот один из вариантов схемы расщепления для системы (1), состоящий из 3-х шагов:
ди' . ди'
— + с Ах — = ^ ,
д* дх
ди" . ди" „ -+ с Ау-= Г
д? ду у
ди" . ди" „
--ь с А2-= Г ,
а & 2
(3)
где Г = Гх + Гу + Г2, Г =
(
0,0,0, -
аЕ„
V
Л
, МВх, -МВХ
2 у
Г =
у
0,0,0,-М Ву,-
°Еу ^ -у, М В
' х
Л
(
Г =
0,0,0, МВ2, -МВг, -
Л
0 J
В качестве начальных условий для каждой системы уравнений из выражения (3) берутся значения, рассчитанные в результате предыдущего шага расщепления. Сохранить второй порядок аппроксимации по времени в схеме расщепления (3) можно с помощью циклического изменения порядка выполнения шагов расщепления, например, выполняя сначала в следующем порядке шаги по пространственным направлениям: х у 2; у х 2 ; 2 х у; х 2 у; у 2 х; 2 у х. Обоснование этого утверждения содержится, например, в монографиях [8, 9].
Рассмотрим первую систему из (3), которая интегрируется на шаге
расщепления по направлению х. Представим матрицу А х в виде -1
Ах = Л, где Лх — диагональная матрица, на диагонали которой стоят
собственные числа матрицы Ах , — матрица, столбцы которой есть правые
собственные векторы матрицы А х , определенные с точностью до множителя, а -1 -1 0 есть матрица, обратная к . Строки матрицы 0 есть левые собственные
векторы матрицы Ах . Эти матрицы можно взять в виде:
Л х =
(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
>вх =
(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
Введем вектор-столбец характеристических переменных w = 0 и,
компоненты которого заданы формулами щ = Вх, w2 = В у - Ег, щ = В 2 - Еу,
щ = Ех , щ = В 2 + Еу, щ = В у + Е . Умножая первую систему из (3) на матрицу -1
0 слева, получим систему из независимых скалярных уравнений для характеристических переменных:
ч
/
т
\
\
дж1 дг = о, дж4 дг
дж 2 + с дж 2
д г дх
дж3 д 1 - с дж3 дх
— Ж , (4)
££0
дж5 дж 5
>ж' 1Г+с аТ=Мж (5)
дж
/ж , —- с —^ = -Мж . (6)
21' & дх у 1 ()
Отметим, что в силу первого уравнения из (4) на шаге расщепления по х постоянны ж и правые части уравнений (5) и (6) и что в силу второго уравнения из (4) для ж4 имеется очевидное аналитическое решение. Таким
образом, численное интегрирование первой системы из (3) сводится к численному интегрированию уравнений (5), которые описывают волны, бегущие слева направо вдоль оси х, и уравнений (6), которые описывают волны, бегущие справа налево.
Для этих уравнений предложено большое число разностных схем [7, 8]. Существует несколько явных монотонных схем, которые имеют 3-й порядок точности по пространству и 1-й порядок точности по времени, например, схема, предложенная в работе [10]. Тестовые расчеты показали, что для получения приемлемого качества численного решения необходимо использовать схемы повышенного порядка точности. Авторами данной работы была предложена следующая явная гибридная схема, которая имеет 2-й порядок точности по времени и 3-й по пространству.
Рассмотрим построение этой схемы. Пусть задана равномерная сетка по времени и равномерная пространственная сетка в декартовых координатах, целые и полуцелые узлы которых заданы соотношениями ги = го + п т,
Х= хо + г Их, у,= Уо + ]Ну, гк= + х1+У2 = х, + 0' +1/2) Их,
Уу+1/2 = У, + О + 17 2) 2к+1/2 = 20 + (к + 1/ 2) hz, в которых Их , Иу , И2 — шаги сетки по осям х, у, 2, а т — шаг по времени. Будем использовать стандартные обозначения для значений функции ж в узлах сетки
ж(гп,xi,у,,) = жпрг ,к . Явные разностные схемы для первого из уравнений (5)
записываются в виде:
п+1 п
- Ж2,г,= -СХ',,,к (+1/2,,,к - ж1г-1/2,,,к ) + тМу,г,, (7)
где С г у к = С} Т / Их > 0 есть число Куранта в направлении х в узле сетки
с индексами г',,,к. Порядок точности схемы (7) и ее монотонность
обеспечиваются способом вычисления значений функции ж2 в полуцелых по х
узлах сетки. Для записи этого способа обычно используют параметр, называемый анализатором гладкости функции в направлении х и заданный формулой
Ш = , (8)
П П ' 4 7
^2,¿+1,у ,к - ,] ,к
и функцию Ш), называемую ограничителем потока. Формулы для
п
вычисления ¿+1/2 . к записывают в виде:
п п
^+1/2, Лк = Лк + (1 - Сх4ш „ )в(Щ "" ^ . (9)
Для схемы Лакса — Вендроффа С(Ш) = 1, для схемы Бима — Уорминга О(Ш) = Ш. Эти схемы имеют 2-й порядок точности по времени и по пространству. При О(ЩШ) = 0 получается монотонная схема 1-го порядка точности. Как показано в работах [11, 12], для монотонности схемы (7), (8) необходимо выполнение следующего условия на ограничитель потока:
0 < а(Ш) < шт(2, Ш+1Ш |). (10)
Из формулы (10) следует, что схема Лакса — Вендроффа монотонна при выполнении условия Ш>0.5, схема Бима — Уорминга монотонна при выполнении условия 0 < Ш < 2. Эти схемы являются немонотонными, однако области монотонности этих схем перекрываются.
В предложенной авторами гибридной схеме осуществляется переключение между шестью схемами, которые имеют 2-й порядок точности по времени и 3-й по пространству и которые монотонны при выполнении условий переключения. Эти условия записываются в следующем виде:
1. Если ]к - ]к = 0, то полагаем О(Ш) = 0 .
1 - С
2. Если н?, м . к - н?,, . к * 0 и Ш< 0, то полагаем О(Ш) = х''' ]'к
2
3. Если н?, г+1, j к к - н?, г, л к * 0 и Ш>0, то
а) при 0 <Ш< (1 - Сх^к )/(3 - С^ ) полагаем С(Ш) = Ш(2 - Сх^к ),
б) при (1 - Сх, , ,, к ) / (3 - Сх, , ,, к ) < Ш < (2 - Сх, , ,, к ) / (5 - Сх, , ,, к ) полагаем О(Ш) = Ш(1 + Сх,,лк ) /2 + (1 - Сх,,лк ) / 2,
в) при (2 - Сх,,м ) / (5 - Сх,,м ) < Ш < (4 + Сх, ,^ ) / (1 + Сх,,м ) полагаем О(Ш) = Ш(1 + Сх,, jЛ ) /3 + (2 - Сх,, jЛ ) / 3 ,
г) при (4 + Схил )/(1 + Сх,,Лк ) <Ш полагаем О(Ш) = (3 + )/2.
Отметим, что пункты «а» и «б» являются уточнениями схемы Бима — Уорминга, а пункт «г» является уточнением схемы Лакса — Вендроффа до 3-го порядка точности по пространству. Пункт «в» является схемой 3-го порядка точности по пространству, использующей центральные разности. Прямой
проверкой можно убедиться, что условие монотонности (10) выполняется для пунктов «а», «б», «в» и «г». Схема для второго уравнения в (5) полностью аналогична изложенной выше.
Рассмотрим схему для первого из уравнений (6). Она записываются в виде:
- м". ■1г= С ■ Лж" /9 - .,)-тМ .М -г, (11)
3,1 ,},к 3,1,},.к х,1,},к\ 3,1+\12,},к 3,1-112,},к / х,1,},к 1,1,},к ' V /
который аналогичен (7). Значения функции ж3 в полуцелых по X узлах сетки вычисляются по формуле:
п п
М3",1+1/2,},к = М3й,1+1,],к + (1 - С] )0(Ш) М3,1М3,1,},к , (12)
в которой анализатор гладкости функции ж3 в направлении X задан формулой
п п
^ _ м3,1+2,},к - М3,1+1,},к
П П ^ '
,1+1,},к - М ,1,},к
Условия переключения для ограничителя потока функции С(^) точно такие же, как в изложенной выше схеме для первого уравнения в (5). Схема для второго уравнения в (6) полностью аналогична изложенной схеме для первого уравнения в (6).
Послойный переход для второго уравнения в (4) осуществляется по формуле:
<!},к = ехр(-а;., },к / ад, },к) м4,1,} ,к,
соответствующей аналитическому решению этого уравнения.
Послойные переходы для второй системы из (3), которая интегрируется на шаге расщепления по направлению у, и для третьей системы из (3), которая интегрируется на шаге расщепления по направлению 2, выполняются аналогично описанному переходу на шаге расщепления по направлению X .
Заключение
Авторами были проведены многочисленные тестовые расчеты, в которых численные решения задач о распространении сигнала от точечного дипольного излучателя в однородном проводнике и в однородном диэлектрике, полученные с помощью представленной в данной работе схемы, сравнивались с точными аналитическими решениями этих задач. Оказалось, что численное решение очень хорошо воспроизводит форму и амплитуду сигнала. На расстоянии от излучателя в пределах 4-5 длин волн амплитудная ошибка численных решений не превышала 2 %. Фазовая ошибка также оказалась очень мала. Численные решения одинаково хорошо воспроизводили как сигналы, состоящие из одного импульса, так и гармонически колеблющийся сигнал. При использовании разностных схем, имеющих 2-й порядок точности по времени и пространству, амплитудная ошибка численного решения увеличивалась до 5 %. Таким образом, предложенная в данной работе схема позволяет получать
численные решения уравнений Максвелла, обладающие более высокой
точностью, чем решения, полученные с помощью схем с порядком точности
меньше 3-го по пространству.
Литература
1. 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. 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. DOI: 10.1007/s10712-009-9063-5.
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. DOI: 10.1109/TAP.2007.897138.
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. DOI: 10.1109/TAP.2011.2167940.
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. DOI: 10.1109/TAP.2009.2037706.
6. Семенов А. Н., Смирнов А. П. Численное моделирование уравнений Максвелла с дисперсными материалами // Математическое моделирование. 2013. Т. 25, № 12. С. 19-32.
7. Куликовский А. Г., Погорелов Н. В., Семенов А. Ю. Математические вопросы численного решения гиперболических систем уравнений. М.: Физматлит, 2001. 608 с.
8. Бисикало Д. В., Жилкин А. Г., Боярчук А. А. Газодинамика тесных двойных звезд. М.: Физматлит, 2013. 632 с.
9. Ковеня В. М., Яненко Н. Н. Метод расщепления в задачах газовой динамики. Новосибирск: Наука, 1981.
10.Вязников К. В., Тишкин В. Ф., Фаворский А. П. Построение монотонных разностных схем повышенного порядка аппроксимации для систем уравнений гиперболического типа // Математическое моделирование. 1989. Т. 1, № 5, С. 95-120.
11.Harten A. High resolution schemes for hyperbolic conservation laws // J. Comp. Phys. 1983. W. 49. P. 357.
12.Sweby P. K. High Resolution Schemes Using Flux Limiters for Hyperbolic Conservation Laws // SIAM J. Numer. Anal. 1984. W. 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
