Локально-одномерная разностная схема для электродинамических задач с заданным волновым фронтом
(Locally one-dimensional finite-difference scheme for the electrodynamic problems with given wavefront
Preprint, Inst. Appl. Math., the Russian Academy of Science)

ВВЕДЕНИЕ

Математические модели электромагнитных полей, генерируемых при ионизации больших объемов газа импульсными потоками фотонов, могут основываться на трехмерных уравнениях Максвелла с начальными условиями на фронте ионизации[1]. Примером такого процесса является образование электромагнитного поля тормозным гамма-излучением, образующимся при воздействии электронов ускорителя на мишень. Если длительность импульса фотонов мала по сравнению с характерным линейным размером заполненного газом объема, то плотность стороннего тока комптоновских электронов и вторичная ионизация образуются в узком слое за фронтом гамма-излучения [2]. Этот факт и обуславливает необходимость выделения переднего фронта. Мишень является, по сути, изотропным источником фотонов. Поскольку длины пробегов квантов существенно превышают размер мишени, существенные участки фронта их потока имеют форму, близкую к сферической. Электромагнитное поле также имеет сферический фронт.

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

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

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

1.Дифференциальные уравнения

Рассмотрим уравнения Максвелла в системе единиц CGS:

,                 ,

где – пространственные координаты,  – лабораторное время, ,  – относительная диэлектрическая и магнитная проницаемость среды соответственно,  – скорость света в вакууме.

Функция  такова, что неравенство  определяет область за фронтом электромагнитного поля. Правая часть уравнений есть плотность электрического тока , где – проводимость среды, а  – плотность стороннего тока.

Можно показать, что в указанных переменных для величины:

из уравнений Максвелла следует соотношение:

                                                                       (En)

Если рассматривать величину  в качестве плотности «энергии» электромагнитного поля, то уравнение (En) выражает закон ее сохранения. Физическая плотность энергии строго положительна для любого отличного от нуля электромагнитного поля. Это выражается в том, что для задачи Коши, когда формально , . В рассматриваемом случае величина, поскольку может обращаться в ноль для ненулевых значений  и  при  .

В декартовой системе координат  с началом координат в центре источника возмущения уравнения Максвелла имеют вид:

где  .

2.Дифференциально-разностная схема

Пусть– область изменения пространственных переменных, в которой рассматривается решение уравнений Максвелла. Введем разностную сетку по переменной:

;

;

Разностную сетку для переменных  введем аналогично.

Выберем параметры сетки так, чтобы разрывы коэффициентов системы уравнений Максвелла разместились так, как показано в таблице 1.

Таблица 1 – Размещение разрывов коэффициентов

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

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

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

;

;

.

Компоненты напряженности магнитного поля  разместим в центрах граней параллелепипедов.

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

   (Hx)

где:

,                                                (Bx)

  для  , если  считать,  что.

Там, где это не вызывает недоразумений, название пространственной переменной в обозначениях для шагов сетки и значения пространственных индексов у значений сеточных функций будем опускать. Соотношения (Hx), (Bx)  удобно записать в виде:

Легко проверить, что справедливо тождество

.                                                          (cBx)

Аналогично определим усредненные на разрыве значения сеточных функций  и:

       (Hy)

;                                                (By)

        (Hz)

,                                                 (Bz)

для которых справедливы тождества, аналогичные (cBx).

Дифференциально-разностные аналоги построим интегрированием уравнений Максвелла по площадкам, определенным в таблице (2):

Таблица 2

Интегрирование проведем, используя теорему о среднем и квадратурные формулы. Квадратурные формулы будут определены самим видом дифференциально-разностных уравнений.

1. Дифференциально-разностный аналог уравнения:

где: ,

.

Здесь и далее используются следующие символы.

Символ  обозначает следующее средневзвешенное значение сеточной функции, заданной в узлах сетки с дробными индексами  в точке сетки с индексами 

.

Символ  обозначает следующее средневзвешенное значение сеточной функции, заданной в узлах сетки с дробными индексами  в точке сетки с индексами 

2. Дифференциально-разностный аналог уравнения

где:  

 .

3. Дифференциально-разностный аналог уравнения:

где: ;

4. Дифференциально-разностный аналог уравнения:

где: ,

5. Дифференциально-разностный аналог уравнения

где ;

.

6. Дифференциально-разностный аналог уравнения:

где: ;

.

Далее будет показано, что для обеспечения положительной определенности разностного аналога энергии в дифференциально-разностных уравнениях –, вместо  следует написать

вместо 

а вместо

3.Дифференциально-разностная теорема об изменении энергии

Умножим и просуммируем дифференциально-разностные уравнения результаты умножений в соответствии с таблицей 2.

Таблица 2

Получим дифференциально-разностный аналог уравнения скорости изменения «энергии» электромагнитного поля:

,

где - разностный аналог мощности тепловых потерь, мощности сторонних токов, а  потока энергии через поверхность расчетной области 

Из этого выражения очевидно, что

.

С учетом определения  справедливы следующие равенства:

и  для  :

где . 

где

где

Можно показать, учитывая эти равенства и определения величин   , что

, 

где

Здесь и далее   ,  ,  .

Справедливы следующие тождества:

,              ,

где ,    ;

,             ,

где ,    ;

,              ,

где ,    ;

, ,

, ,

где  ,

,;

, ,

,  ,

где ,

, ;    

, , 

, , 

где  ,

,

.

Тогда

,

,

,  

,

,

,

и

,

,

.

Если учесть тождество:

,

то

Здесь и далее   .  Аналогично,

и

На основании вышеизложенного и после ряда преобразований, аналогичных преобразованиям для получения  (En),  выражение для плотности разностной энергии элементарного объема  приводится к виду:

,

где:

очевидно, является разностным аналогом плотности «энергии» ЭМП, так как   и  аппроксимирует ее со вторым порядком, а  является неотрицательной величиной второго порядка размера элементарного объема,  что сразу следует из приведенного ниже выражения 

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

3.Решение сеточных уравнений

Запишем исходную систему уравнений Максвелла в матричном виде:

        (3.1)

где u – вектор значений напряженности электрического и магнитного поля

;;;;

;;

Введем неравномерную сетку по переменной  с узлами. Обозначим        τ_n = t_n+1 – t_n. Разобьем интервал [t_n , t_n+1] на шесть подинтервалов одинаковой длины τ_n /6. Будем считать матрицы C и F на интервале [t_n , t_n+1] постоянными. На каждом из подинтервалов заменим систему дифференциальные уравнения своей системой разностных уравнений. Обозначим .

Введем вспомогательную индексацию временных слоев (рисунок 3.1):

₀

t₁                t₂                t₃                t₄                t₅                t₆     

t_n                                              t_n+1/2                                    t_n+1

Рисунок 3.1 – Индексация промежуточных временных слоев.

В соответствии с введенной индексацией обозначим

u⁰ = u(t₀);u¹ = u(t₁),…, u⁶ = u(t₆),  u⁰¹ = (u⁰+u¹)/2,u¹² = (u¹+u²)/2,…      ,и т.д.

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

                                       (3.1.1)

                                       (3.1.2)

                                       (3.1.3)

                                        (3.1.4)

                                        (3.1.5)

                                        (3.1.6)

Система уравнений (3.1.1-6) является симметричной разностной схемой. Очевидно, что ни одно из уравнений схемы не аппроксимирует уравнение (3.1). Сложим уравнения (3.1.1-6):

Разлагая все функции, входящие в дифференциальную постановку, в ряд Тейлора по времени относительно t= t₃ = t_n+1/2, нетрудно убедиться, что данное уравнение аппроксимирует  дифференциальное уравнение (3.1) с порядком O(τ²) при постоянном шаге по переменной.  

Рассмотрим более подробно уравнение (3.1.1). Оно представляет собой систему шести разностных уравнений относительно компонент вектора u:

                                                                                (3.1.1.1)

                                         (3.1.1.2)

                               (3.1.1.3)

                                                                                (3.1.1.4)

                                   (3.1.1.5)

                                  (3.1.1.6)

Уравнения (3.1.1.1) и (3.1.1.4) являются явными.

Система уравнений (3.1.1.2), (3.1.1.6) при каждом j и   k  распадается на подсистемы, зависящие только от значений функций при данных  j и k с неизвестными E_y¹ и H_z¹ на верхнем слое. Эта система решается прогонкой по переменной x.

В уравнения (3.1.1.3) и (3.1.1.5) входят неизвестные E_z¹ и H_y¹ с верхнего временного слоя. Она также решается прогонкой по x. Для остальных подинтервалов аналогично.

4.Тестирование разностной схемы.

Эффективность разностной схемы и метода решения разностных уравнений проверена на следующей серии тестов.

1. Идеальный проводник на внешней границе.

Расчетная область:, ,.

Плотность тока:

,  ,.

Проводимость , диэлектрическая проницаемость.

Точное решение:

, , ,

,,

2. «Финитный» тест (цилиндрическая симметрия). Рассмотрим аксиально-симметричную задачу в цилиндрической системе координат (ρ,φ, z). Зададим электрическое поле:

,  ,  где  ,

Здесь T(t) – временная зависимость, а Ψ_ρ и Ψ_z  – некоторые финитные функции.

Тогда магнитное поле.

Такое электромагнитное поле получается при следующей плотности тока:

В качестве финитных функций возьмем следующие (на примере Ψ_z ):

,

где

Данный тест составлен В.И. Турчаниновым.

3. Тест с постоянной проводимостью. Рассмотрим аксиально-симметричную задачу в сферической системе координат (ρ, θ, φ). Пусть плотность тока имеет только радиальную компоненту.

Аналитическое решение имеет вид:

где ;

В формулах для электрического поля использованы следующие обозначения:

;

; 

Тест составлен С.В. Петропавловским.

4.Электрический диполь. Рассмотрим аксиально-симметричную задачу в сферической системе координат (r, θ, φ) с плотностью стороннего тока:

,   

Аналитическое решение имеет вид:

,

,     

где , , .

Заключение

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

ЛИТЕРАТУРА

1. Турчанинов В.И. Численная методика решения трехмерных уравнений Максвелла в сферических переменных в неоднородной диссипативной среде с выделением  переднего фронта. М. Препринт ИПМ им. М.В. Келдыша РАН, 1993 №18

2. Гайтлер В. Квантовая теория излучения. – М.: ИЛ, 1956.