CC BY
27
УДК 519.5, 537.8
УРАВНЕНИЯ МАКСВЕЛЛА В ЗАДАЧЕ РАССЕЯНИЯ СВЕТА И ЧИСЛЕННЫЕ МЕТОДЫ ИХ РЕШЕНИЯ
© 2010 г. И.Э. Наац1, Э.И. Наац2
1Ставропольский государственный университет, 1Stavropol State University,
ул Пушкина, 1, г. Ставрополь, 355009, Pushkin St., 1, Stavropol, 355009,
info@stavsu.ru info@stavsu.ru
2Институт оптики атмосферы СО РАН, 2Institute of Atmospheric Optics SB RAS,
пл. Академика Зуева 1, г. Томск, 634021, Academic Zuev Sq., 1, Tomsk, 634021,
contact@iao.ru contact@iao.ru
Рассмотрено применение метода взвешенных невязок к численному решению уравнений Максвелла. Получены общие формулы метода в случае произвольных базисных и поверочных функций. Рассмотрена реализация алгоритма с применением конечно-элементной аппроксимации. Работа алгоритма проиллюстрирована на примере задачи рассеяния плоской электромагнитной волны на проводящей диэлектрической сфере.
Ключевые слова: уравнения Максвелла, метод взвешенной невязки, вычислительный алгоритм, вычислительный эксперимент, задача рассеяния плоской электромагнитной волны.
Application of the method of weighted residuals for the numerical solution of Maxwell's equations is considered. We obtain general formulas of the method in the case of arbitrary base and calibration functions. We consider the implementation of the algorithm using the finite element approximation. The algorithm is illustrated by the example of the problem of electromagnetic dispersion by a conducting dielectric sphere.
Keywords: Maxwell's equations, weighted residuals method, computational algorithm, computing experiment, the problem of dispersion plane electromagnetic wave.
Задачи рассеяния света несферическими частицами представляют собой одну из интенсивно развивающихся областей физики. К ним относится изучение рассеяния света кристаллическими облаками в атмосфере. Особое значение при этом уделяется изучению рассеяния излучения облаками посредством наземных и космических лидаров (один из масштабных проектов - CALIPSO [1]). Среди задач, стоящих перед исследователями, особое место занимает разработка корректных моделей, которые необходимы для оценок вклада кристаллических облаков в рассеяние излучения лидаров [2]. Модели строятся на основе знаний о рассеивающих свойствах частиц, составляющих облака, причем требования к точности рассчитываемых характеристик постоянно возрастают. Проблемы в основном определяются размерами рассеивающих частиц по отношению к длине волны падающего света. Чем больше размеры частиц, тем больше требований возникает к вычислительным ресурсам и, следовательно, к алгоритмам численного решения задачи рассеяния.
В настоящее время имеется много подходов к численному решению таких задач [3]. В случае больших размеров частиц распространение получили методы геометрической оптики [4]. Однако точность таких методов недостаточна, особенно для получения значений такой характеристики рассеяния, как матрица Мюллера. Теоретически наиболее точное решение может быть получено непосредственно численным решением уравнений Максвелла. Из наиболее известных и хорошо исследованных методов численного решения уравнений Максвелла является метод FDTD (finite difference time domain) [5]. Однако и у этого метода есть свои ограничения, в основном связанные с требованиями к вычислительным ресурсам, особенно при расчетах для частиц с большими размерами по отношению к длине волны.
В данной работе авторами предлагается численный метод решения уравнений Максвелла в задаче рассеяния электромагнитных волн в рамках вариационного подхода (метод взвешенных невязок [6]). Построенный на его основе алгоритм численного решения обладает рядом важных свойств, которые делают его эффективным в программной реализации, прежде всего в случае больших размеров рассеивающих частиц.
Постановка задачи рассеяния света на основе уравнений Максвелла
Координатная форма уравнений Максвелла, связывающих рассеянное поле (E, H) с падающим полем (E ^, H ^), записывается так [7]:
дЕ дЕ„
jko ^ S'EX =
jk0 ^ е'Еу =
sy
S^S s
SyS
s
д~ дн
j^JUL S'EZ =
dy dz
dH x dH z
dz dx
dH y dH x
y - jk0 (e- iE>;
- jk0 (e'- i)EJ);
дх dy
ßo SJSLH x =-Ц- ^;
dy dz
- jk0 (e- 1)e
(i) (2)
(3)
(4)
jk0^Hv =--^ + ■
dz
dx
jko ^tf. = -
s^ ~ dx
dEy+ дЕх
dy
В уравнениях (1) - (6)
CTV .
jrne
H = щ H ; щ = ™
, e a e =--j-
(5)
(6)
(7)
(8) (9)
ax (x)=ljl ax,max , Kx
где е0 и - электрическая и магнитная проницаемости среды; а - величина проводимости в пределах зоны РМЬ; к0 - волновой вектор. Для параметров в задаче (7) приняты следующие рекомендации [5]:
( Лт
= 1 + ^,тах• (I0)
Внутри зоны РМЬ величина ах меняется от нуля до сгх тх , вне - равна нулю. Соответственно, кх внутри
зоны РМЬ возрастает от единицы до кхтах , вне - равна
единице. Параметр т в (10) выбирается в диапазоне 3<т<4. Остальные параметры подбираются в ходе вычислительного эксперимента с учетом эффектов погашения приходящих и отраженных волн. Для поля падающей плоской волны в расчетах использовались рекомендации из [7]. Считалось при этом, что падающая волна может быть представлена в сферических
(i) = (Ев-0 + Еф • ф)ехр((г '• P)/c + R/c),
координатах: Е{) = \Ев
Н{) = {еф-§- Ев- ф)ехр((г ' • г)/ с + И./с) и распространяется в направлениях, определяемых сферическими координатами (в, ф). Единичный вектор г направлен из начала координат в направлении (в,ф),
вектор г направлен из начала координат в некоторую точку внутри области расчета; Я - произвольное расстояние. В декартовой системе координат амплитуды падающей волны представляются в виде
Е() = Ев с08(в)с08(ф) - Еф 8т(ф),
ЕУ) = Ев с08-в)8т(ф) + Еф соб(ф) , Е^ = -Ев 8т(в),
Н^ = Ев sin(ф)+ Еф а^(в)о^(ф), НУ) = -Ев cos(ф) + Еф о^(в>т(ф), Н= -Ев sin(ф).
Для представленной выше математической модели процесса рассеяния электромагнитных волн на частицах произвольной формы авторами разработан в рамках вариационного подхода численный метод (алгоритм), создано программное обеспечение и проведен вычислительный эксперимент. Соответствующие результаты приводятся ниже.
Численное решение методом взвешенных невязок
Алгоритм метода взвешенных невязок применительно к решению рассматриваемой задаче можно описать следующим образом [6]. Пусть имеется некоторый линейный оператор Ь (в частности дифференцируемый) и
s
a
и
y
z
sx =Kx +
e
0
s
x
f — функция «источника». Тогда широкий класс задач может быть представлен в виде уравнения
Lu = f . (11)
Невязкой уравнения (11) принято считать величину г = Lu - /. (12)
Пробное решение un по методу взвешенных невязок подбирается таким образом, чтобы выполнялись следующие условия для невязки (12):
|(¿ип - Г)= о. (13)
□
Здесь О представляет собой область, в которой существует решение уравнения (11), а функции называются поверочными функциями. Далее пробное решение может быть представлено как разложение по некоторым базисным функциям Ф^ в виде
un = Z ckФk . k=1
(14)
Базисные функции подбираются таким образом, чтобы удовлетворить граничным условиям точно или приближенно [6]. Подставляя (14) в (13) и используя свойство линейности оператора L, получим
п
2 ск | (¿Ф к Г .
к=1 □ □
Очень часто форма оператора L позволяет провести интегрирование по частям, что приводит к снижению требований к гладкости базисных функций [6]. В итоге мы приходим к системе линейных уравнений для определения неизвестных коэффициентов
п
2 а1кск , ^ = |fT/d□, ал = |(¿Фк.
к=1 □ □
Применение метода взвешенных невязок к системе (1) - (6) проиллюстрируем на примере первого уравнения (1). Введем пробные решения:
Ei"](x, y, z ^¿c^ )ф k (x, y, z),
к=1
Hyn)(x, y, z blc^ к (x, y, z), к=1
H zn)(x, y, z )=l4Azb к (x, y, z) .
(15)
(16) (17)
к=1
„(ет) c(hy) c (hz)
Неизвестные коэффициенты с( , ск , ск
комплексные, а базисные функции - вещественные, удовлетворяющие однородным граничным условиям на границе области (обращающимся в нуль на границе). Далее введем обозначение для функции источника: /х = А 0?'-1)4г).
Представим функцию источника в виде разложения по тем же базисным функциям
fLn)(x,y,z)=iZc^к(x,y,z) .
(18)
к=1
Коэффициенты с(ех) в (18) - комплексные, могут
быть найдены путем решения задачи аппроксимации методом наименьших квадратов. Однако, как будет показано ниже, при выборе локальных базисных функций в этом нет необходимости. Далее подставим (15) - (18) в уравнение (1) и, применяя условие (13), получим
ZcS«)
к=1
А i
SySz „i
е'Ф к T{dü
ü s x
к=1
- ZcH |Фк^idü]-iZc^y)i - J-^Фкdü I. (19)
к=1
=icih"'[-J^^L Ф Idüj-
-y
к=1
dz
Второе и третье слагаемые в правой части получаются интегрированием по частям и учетом граничных условий. При этом подразумевается, что область интегрирования О представляет собой куб или параллелепипед, содержащий рассеивающий объем.
Аналогичным образом вводятся пробные решения для остальных компонент поля и подставляются в остальные уравнения (2)-(6) с последующим вычислением интегралов. Приведем окончательные выражения для элементов матриц и получающихся алгебраических уравнений для определения неизвестных коэффициентов. Введем следующие обозначения для элементов матриц:
ОТ) = А J—е'Ф к Tidü; flW = ]кй J^ Ф к Tidü; (20)
ü sx ü Sx
а,к = ih J—е'ФкTdü; a(hy) = jh0 J-^Фк^dü ; (21)
ü Sy
ü Sy
ак) = A J—е'Ф к %dü; a^ = j^ J^ Ф к ^dü ; (22)
ü sz ü Sz
bix)=-Jd^LФkdü; ьк)=-Jd^^Фkdü ; (23)
dx
—
ü
i \ -T
Ьк) = - J-T^Ф^ ; dlk = |ФкTidü. (24)
□ бг □
Система алгебраических уравнений (19) в матрично-векторных обозначениях будет иметь следующий вид:
A(ex )c(ex) = B(у)c(hz) - B(z )c(^y)- Dc(,ex)
A(ey )c(ey) = B(z)c(hx)- B(x)c(hz) - Dc(ey)
A(ez )c(ez) = B(x)c(^y)- B(y)c(hx) - Dc(iez)
A(hx)c(hx) = -B(y)c(ez + B(z )c(ey)
= -B(z )c(ex) + B(x)c(ez)
A(hz)c(hz) = -B(x)c(ey) + B (y)c(ex)
Более компактно ее можно представить с помощью вспомогательных матриц
А« =
АН 0 0 A(ey)
0 0 A(ez)
0 B(z) -B(y)
B = - B (z) 0 B (x)
B (y) -B(x) 0
A(h) =
A(hx) 0 0 A(hy)
0
0
0 0
A(hz)
А =
A(e) B - B A(h)
следующим
образом: АС = Р, где введены обозначения: С = (с(ет) с(еу) с(ег) с(Ах) с(Ау) с(Аг)) Р = (рс Ы, Пс('еу), Рс(е) ,0,0,0).
Как показали численные эксперименты, наиболее эффективным методом (по крайней мере, для выбранных базисных функций, описанных ниже) может быть неявная итерационная схема:
С(п+1) - С(п)
■ = F - AC
(и+1)
или
ü
ü
ü
n
т
(Е + т • А)С{и+1)= С{и) + т Р, (25)
где Е - единичная матрица; т - вещественный итерационный параметр. Начальное приближение для коэффициентов выбирается в соответствии с нулевыми условиями. Итерационный параметр подбирается в ходе вычислительных экспериментов. Программная реализация алгоритма (25) может быть выполнена с применением методов распараллеливания алгоритмов.
Теперь необходимо задать конкретный вид базисных и поверочных функций. Будем выбирать базисные и поверочные функции в виде суперпозиции одномерных функций, например:
ф ^ г) = ф (х) • (у) • фк (г) • (26)
Если область решения представляет собой куб, и количество базисных функций вдоль каждой оси одинаково и равно М, тогда для индекса ^ в (26) справедлива
формула: £ = (к-1)-М2 + (/-1) М +1 , 1 <£< N, N = Мъ.
При таком выборе вида базисных и поверочных функций вычисление интегралов (20)-(24) упрощается. Как показали вычислительные эксперименты, для поверочных функций можно выбирать особую форму, которая обеспечивает большую устойчивостью численного решения.
Разобьем область решения вдоль каждой пространственной оси на интервалы в общем случае неравной длины (неравномерная сетка). В расчетах были выбраны локальные одномерные базисные функции вида
V (х) =
0; х < хк_j; х > х
х - х
k-1
х1- х
k-1
хк+1 х
хк+1 хк
k+1
хк-1 < х < хк
хк < х < хк+1 •
Для поверочных функций форма выбиралась следующим образом [6]:
¥к{х)=фк{x)±a•yk{х), 0<а < 1, (27)
7 к (х) =
0; х< хк-1; х> хк+1
(х - хк -1 Хч - х).
хь х.
хк-1 < х < хк
(28)
к-1
(хк+1 - х Xх - хк ) .
новое число. Методика расчета индикатрисы рассеяния Р(в) полностью соответствует [5]. Область расчета представляла собой куб с длиной ребра, равной 3, окружающий рассеивающую проводящую сферу единичного радиуса с показателем преломления п = 1,311+/-10-4 (типичное значение для ледяных кристаллов). Падающая электромагнитная плоская волна распространяется в
направлении оси Ог слева направо с линейной поляризацией вдоль оси Ох. Число базисных функций бралось одинаковым по всем трем измерениям с постоянным интервалом разбиения. Во всех расчетах параметр а = 0,5. Для расчетов с параметром ка = 5 число локальных одномерных базисных функций равнялось 71, а для параметра ка = 24 составляло 151. Соответственно, выбирались параметры сттах = 10 и ктх = 5 для первого расчета и сттах = 50 и ктах = 5 для второго. Ширина
зоны РМЬ составляла 10 интервалов для обоих расчетов.
На рис. 1, 2 приведены результаты расчетов для значения параметра ka=5, а на рис. 3, 4 - для значения параметра ка = 24.
Рис. 1. Индикатрисы рассеяния, полученные по теории Ми и алгоритма ^ОТБ для значения параметра ка = 5
хк+1 хк
Выражения (27) и (28) предполагают, что используется зеркально-симметричная схема относительно начала координат, расположенная в центре куба (параллелепипеда), окружающего расчетную область. Тогда, если поверочные функции расположены левее начала координат, то используется «+», иначе - «-».
Пример численного решения
Для иллюстрации работы описанного выше алгоритма численного решения уравнений Максвелла рассмотрим задачу рассеяния плоской электромагнитной волны на проводящей сфере для двух значений параметров: ка = 5 и ка = 24, где а - радиус сферы; к - вол-
-1,5 "-1,5
-1,0 -0,5
0,0 0,5
Z
Рис. 2. Рассчитанное распределение амплитуды полного электрического поля Ех в плоскости уОг для значения параметра ка = 5
1,5
1,0
0,5
хк < х < хк+1
Y 0,0
0,5
1,0
10*
10-" -I-1-1-1-1-1-1-1-1-
О 20 40 60 80 100 120 140 160 180 в
Рис. 3. Индикатрисы рассеяния, полученные по теории Ми и алгоритма ^ОТБ для значения параметра ка = 24
Z
Рис. 4. Рассчитанное распределение амплитуды полного электрического поля Ex в плоскости yOz для значения параметра ka = 24
На рис. 1, 3 представлены индикатрисы рассеяния, рассчитанные по теории Ми и разработанному алгоритму WRFD (weighted residuals frequency domain).
На рис. 2, 4 приведены рассчитанные распределения компоненты полного электрического поля Ex в плоскости yOz. Приведенные результаты расчетов наглядно демонстрируют хорошее согласие с теорией Ми. Следует отметить, что в расчетах не ставилось целью достичь полного согласия двух подходов, но лишь добиться приемлемых результатов при минимально необходимых ресурсах компьютера и времени расчета. Приведенные на рис. 1 - 4 результаты расчетов демонстрируют то, как должно возрастать для увеличения точности получаемых результатов число базисных функций при усложнении формы индикатрисы. Обращает на себя внимание усложнение картины распределения поля внутри и около сферы при разных длинах волн падающего излучения.
В целом оба рассмотренных случая позволяют оценить требования к вычислительным ресурсам при уменьшении длины волны и, соответственно, росте значения параметра ka. Однако эти ресурсы заметно меньше, чем в аналогичных конечно-разностных методах, например, FDTD.
Литература
1. The CALIPSO mission: Spaceborne lidar for observation of aerosols and clouds / D.M. Winker [et al.] // Proc. SPIE. 2003. Vol. 4893. Р. 1-11.
2. Боровой А.Г., Кустова Н.В. Разработка оптической модели перистых облаков // Оптика атмосферы и океана. Физика атмосферы : материалы XVI Междунар. симп. Томск, 2009. С. 514-516.
3. Mishchenko M.I., Travis L.D., Lacis A.A. Scattering, Absorption, and Emission of Light by Small Particles. Cambridge, 2002. 690 c.
4. Yang P., Liou K.N. Geometric-optics-integral-equation method for light scattering by nonspherical ice crystals // Appl. Opt. 1996. Vol. 35. P. 6568-6584.
5. Taflove A., Hagness S.C. Computational Electrodynamics: the Finite-Difference Time-Domain Method. Norwood, 2000. 852 с.
6. Флетчер К. Численные методы на основе метода Га-леркина : пер. с англ. М., 1988. 352 с.
7. Kunz K.S., Luebbers R.J. The Finite-Difference TimeDomain Method for Electromagnetics. Boca Raton, FL, 1993. 448 р.
Поступила в редакцию_25 января 2010 г.
