Диэлектрические резонаторы: метод интегральных и интегродифференциальных уравнений Текст научной статьи по специальности «Математика»

ФИЗИКА

УДК 621.372.8.001.24

ДИЭЛЕКТРИЧЕСКИЕ РЕЗОНАТОРЫ:

МЕТОД ИНТЕГРАЛЬНЫХ И ИНТЕГРОДИФФЕРЕНЦИАЛЬНЫХ УРАВНЕНИЙ

М.В. Давидович

Саратовский государственный университет, кафедра радиотехники и электродинамики E-mail: DavidovichMV@info.sgu.ru

Рассмотрены объемные интегральные и интегродифференциальные уравнения, описывающие стационарные собственные моды открытых диэлектрических резонаторов, конкретизированные для цилиндрического случая. Предложены алгоритмы их решения. Численно итерационными методами исследованы собственные колебания Н 01s и Н 011 цилиндрического диэлектрического резонатора.

Dielectric Resonators: The Integral and Integrodifferential Equations Methods M.V. Davidovich

The volume integrodifferential and integral equations which are describing the stationary eigenmodes of open dielectric resonators have been considered and concretized for cylindrical case. The numerical algorithms for its solutions are proposed. The H01<?and H011 eigen oscillations of cylindrical dielectric resonator have been numerically investigated using the iteration technique.

Введение

Методам анализа диэлектрических резонаторов (ДР) посвящено большое число работ (см., например, монографию [1] и приведенный в ней список литературы), при этом использовались различные методы: приближенные эвристические (например, метод магнитной стенки), теория возмущений (разложение по малому параметру) [2], метод частичных областей (МЧО), или сшивания, с получением поверхностных интегральных уравнений (ПИУ) [1], метод ПИУ [3], метод объемного интегрального уравнения (ОИУ) [4]. Из них строгими являются методы ПИУ и ОИУ, причем последний является наиболее универсальным, поскольку позволяет анализировать неоднородные ДР произвольной формы. Однако в литературе он достаточно полно не рассмотрен. В публикации [4] с помощью него анализируется мода H01s однородного цилиндрического диэлектрического резонатора (ЦДР), причем в конечном итоге используется теория возмущений. Из-за сложности определения коэффициентов разложения МЧО для ЦДР так и не был строго реализован, а был заменен приближенным подходом [1], не позволяющим определять радиационные добротности колебаний.

В работе рассмотрены ОИУ и объемно-поверхностные ин-тегродифференциальные уравнения (ИДУ) для произвольных анизотропных и неоднородных ДР, конкретизированные и численно исследованные для ЦДР.

© М.В. Давидович, 2007

1. Интегральные и интегродифференциальные уравнения

Открытый диэлектрический резонатор в виде диэлектрического тела объема V с проницаемостью ё(т), ограниченного замкнутой поверхностью S произвольной формы, описывается однородным объемным гипер-сингулярным интегральным уравнением (ИУ)

E (r )= (W • + k2 )х

х J G (k, r - r ')[г£(г')- €]e (r ')dV(1)

V

В уравнении (1) k = ft) / c, G (k, r ) =

= (4n \r |) 1 exp (- jk\r |), проницаемость £

может, вообще говоря, зависеть от частоты и

быть тензором, соответственно чему £ означает единичный тензор. Далее явные зависимости от частоты (волнового числа k) для скалярной функции Грина G, проницаемости

и электрического поля E будем для удобства опускать. Формальное внесение дифференциального оператора W • = grad • div, действующего в (1) на нештрихованные координаты (координаты точки наблюдения), под интеграл приводит к появлению в ядре сильной неинтегрируемой особенности типа

|r r/| — 3

r — r , что делает невозможным применение кусочно-постоянных функций (или квадратурных формул) для непосредственного численного решения. Получающийся интеграл в обычном смысле расходится, и правую часть (1) следует понимать в обобщенном смысле как распределение (обобщенную функцию) [5], а ИУ (1) можно трактовать как псевдоинтегральное. Операторы, соответствующие ИУ (1), называют также псевдодиф-ференциальными [5-7] .

ИУ (1) определяет комплексные волновые числа k (частоты колебаний) соответствующего ДР, при этом зависимость проницаемости от координат позволяет описывать большое разнообразие структур. Модификация ИУ (1) путем введения поверхностных электрических токов и объемного магнитного тока поляризации делает возможным учет металлических экранов и магнитных свойств ДР. В целях компактности изложения эти задачи здесь не рассматриваются. Если проницаемость однородна и может быть выне-

сена за знак интеграла, ИУ (1) существенно упрощается. Вместо ИУ (1) будем использовать эквивалентное ему интегродифференци-альное уравнение (ИДУ) [8]:

Е (г )= | {к2G (г - г ')[£(т')- €]Ё (г ') +

У

+ VG(г - г')У' • [(е(г')- €)Е(г')]}^У' + (2)

+ £[е(г')- £](('(г')Е(г')^'в(г - г')с1Б

Б

Термин ИДУ будем применять к уравнению, в котором на неизвестную функцию действуют как операторы интегрирования, так и операторы дифференцирования, причем ее производные могут находиться как под знаком интеграла, так и вне него. Последний случай сводится к первому формальным введением дельта-особенности в ядро. В [9] и в ряде других работ под ИДУ понимают уравнение типа (1), в котором дифференциальный оператор действует на интегральный оператор (интеграл). ИДУ (2) является сингулярным с интегрируемой особенностью. Его можно получить, используя, например, формулы Стреттона-Чу [9], или же непосредственно применяя векторные интегральные теоремы и перенося в (1) операции дифференцирования с ядра на подынтегральную функцию [8]. Поскольку нормальная компонента электрического поля имеет скачок на поверхности Б, в поверхностном интеграле в формуле (2) следует взять ее внутреннее значение. Внешнее значение, соответствующее вакууму, будем обозначать символом «+», а внутреннее - соответственно как Е = Е-. Наличие скачка приводит к появлению на Б поверхностной плотности связанного заряда д (г )= [е (г)- /]е(-(? ), г е Б. В случае скалярной проницаемости ИДУ (2) может быть преобразовано в объемно-поверхностное ИУ [4, 9]:

Е (г) = £ {к2 в (г - г ')[е(г')- 1]Е (г')-

У

- £ Х (г ')[Е (г '^'е(г ')^'в (г - г ')}йУ'+ (3) + £ [е(г')- 1](г'(г ')Е (г ')^'в (г - г ')ёБ

Б

поскольку тогда в силу соленоидальности вектора е(г )Е (г) имеем

V • [(е(г) - 1)Е (г )] = -V • Е (г) = е~х (г )Е (г )Уг(г).

В этом случае а(г ) = [^(г)-1] Е,-(г ) = Е^(г )-- Е-(г), т.е. поверхностная плотность заряда

является скалярной величиной и выражается через скачок поля. Поверхностные интегралы в (2) и (3) появляются, если тело имеет резкую границу, т.е. функция &(г) на 5 разрывная и скачком уменьшается до единицы. Если же ё{г) гладкая и плавно уменьшается до единицы в некотором внутреннем приповерхностном слое, то Е+(г ) = Е~(г) и поверхностный интеграл не возникает, однако в указанном слое имеется объемная плотность заряда. Нетрудно показать, что в пределе при уменьшении толщины такого слоя до нуля объемный интеграл по нему эквивалентен поверхностному интегралу от получающегося скачка нормальной компоненты поля. Еще один подход к анализу ДР может быть основан на сингулярном ИУ, получаемом из (1) путем выделения особенности в(0, г) в ФГ в(к, г ) = в(0, г ) + Лв(к, г) и соответствующего ей внеинтегрального члена путем интегрирования по бесконечно малому шару, окружающему точку наблюдения. Такое ИУ не содержит поверхностных интегралов.

Далее будем рассматривать ИУ (3), поскольку оно описывает достаточно широкий класс краевых задач для ДР и имеет по сравнению с (1) ядро с пониженной максимальной особенностью типа производной потенциала простого слоя [10]. Если проницаемость является кусочно-постоянной величиной, то объем V (и соответственно интеграл) разбивается на подобласти, внутри которых Ve = 0. При этом нормальная компонента поля приобретает скачки на границах раздела, и необходимо добавить к поверхностному интегралу соответствующие интегралы по двусторонним поверхностям раздела 5 *. В этом случае ядро объемного интегрального оператора имеет слабую особенность типа

|г — г'| 1. При нахождении точки наблюдения

на поверхности особенность в поверхностном интеграле V'G является интегрируемой, если интеграл понимать в смысле главного значения. При этом интеграл от нее по участку, вырезаемому бесконечно малой сферой с центром в точке r , равен нулю. Поэтому к ИУ (3) можно применять численные методы с использованием кусочнопостоянных объемных и поверхностных элементов (или соответствующие квадратурные формулы).

Рассмотрим цилиндрические ДР, для чего используем представление ФГ в цилиндрической системе координат [11]. Один из способов такого представления заключается в замене переменных x = pcos(^),

y = psin(^) в соответствующих декартовых представлениях, например,

G(p,ty,z\p',ty',z') = (4nR) 1 exp(—jkR), (4) R = (p2 + p'2 — 2pp' cos(^ — <p') + (z — z')2).

Можно выполнить указанную замену и в других декартовых представлениях ФГ, описываемых приведенными в [11] соотношениями, например, в (2.8) и (2.14), или же использовать представление (2.15) (далее все аналогичные ссылки будут соответствовать указанной работе). Однако часто более удобно применять непосредственные представления ФГ в цилиндрической системе (формулы (2.17), (2.18), (2.20), (2.21) и (2.22)). Наличие нескольких видов представлений ФГ создает удобства при вычислении интегралов в матричных элементах, поскольку есть возможность выбрать наиболее подходящую формулу, для которой они вычисляются наиболее просто, а также позволяет решить задачу несколькими алгоритмами. При переходе в цилиндрическую систему координат вектор-

потенциал A и вектор E преобразуются по формулам (2.63), при этом ИУ (3) приобретает вид

С (р, г I р', г0 =

Ер(р,ф, г) = к2£в(г - Г)[г(Г)- 1]х

V 2п

х[е(ґ)е08(^-<р') + еАт')ап(9-^/)}^V/- = |с°8 (9-¥)с(р>9>г 1 р',9г')' ■

(9)

-^(ї' )[Е(Г )Уг(г' )]дС(р/ ) } £ [г(г')- 1](г(г' )Е(г' ))Щ=Г) дБ',

(5)

+ £ |£(Г

Б

Е,

(р,9, г) = к2 £ С(г - Г)[г(Г) - 1]х

V

х [е9(ї' )с°8(9- 9) - Ер(ї' Ьт(9 - 9/)}^'

- £ |г- (г')[Е(г')Уг(г')]дСРд^1^' +

£ [&)- 1](г(г' )Е(г' ))дСррд^ дБ',

(6)

+ £

Б

Е (Р,9’ г ) =

= к2£ С(г - г')[г(г') - 1]Ег (г')м' -

V

- £ ^ (г')[Е(г')У'г(г')]дС(дг- г )+ (7) £ [г(г') - 1]Иг')Е(г'))дС(^7 Г ) дБ

+ £ г(г

Б

В связанных ИУ (5)—(7) для компактности записи использованы координаты г и г , которые следует рассматривать в цилиндрической системе (например, как в (4)). Элемент объема имеет вид dV' = р'йр'йф'йг', а вид элемента поверхности йБ' зависит от координаты точки на цилиндре.

В случае, если поле не зависит от координаты ф, указанные ФГ можно упростить, выполнив интегрирование по углу. Пусть, например, рассматриваются азимутальносимметричные Н-колебания изотропного ЦДР радиуса г0 и высоты к. Тогда отлична

от нуля лишь Еф -компонента, и ИУ (5)—(7) приобретает вид [4]

EФ(P, г) =

= к2 £ (е(р', г )- 1)в(р, г I р , г )р йрйг , (8)

БМ

где БМ - меридианное сечение: 0 <р< г0, Ы < к / 2, а ядро имеет представление

ФГ (9) не зависит от ф и представляется также удвоенным интегралом типа (9) по области (0,п). Действительно, сделаем в (9) замену переменных $ = ф' - ф. Поскольку Я и в периодичны по $ с периодом 2п , значение интеграла (9) по области (- ф, 2п - ф) не зависит от значения ф. В [4] предложено вычислять ФГ (9), разлагая экспоненту в (4) в ряд и выражая рекуррентно получающиеся интегралы через эллиптические интегралы первого и второго рода. Но для ФГ (9) есть и другие представления, например, из (2.18) следует

в (р, г I р , г ) =

1 (хр)3\ (хр,)ехр(-Л/х2 - к21г - г')

4%2 - к 2

а из (2.22) соответственно получим С (р, г I р , г' ) =

(10)

(?(уік2 - у2 А1 (у[к

2

Г р)

4 £ 1^1 (л/ к2 - Г2р']н12)(л1 к2 -ГР х ехр (- ]у(г - г'))йу.

х

(11)

В (11) верхнее значение в фигурной скобке надо брать при р < р', а нижнее -при р > р'. Представления ФГ (9)-(11) удобны для анализа азимутально-симметричных Н0п8 и Е0п8 типов колебаний. Для гибридных НЕтп8 и ЕНтп8 типов ФГ имеют вид (2.18) и (2.22), в суммах которых необходимо оставить только один азимутальный член

р(-]т(ф-ф')) ■

ех

2. Понижение особенности

методом непосредственного интегрирования

Рассмотрим другой метод получения ИУ с пониженной особенностью, заключающийся в интегрировании ИУ (1) по координатам точки наблюдения. Для простоты будем считать проницаемость скалярной величиной и константой. Согласно теореме Гельмгольца любой вектор однозначно мо-

2

жет быть представлен своей соленоидальной и потенциальной частями [12, 13], поэтому разобьем электрическое поле на соленои-дальную и потенциальную части: Е=Ё!+Ёр,

или

Е(г ) = Ух С-УФ .

(12)

Далее индекс я означает соленоидаль-ную, а р - потенциальную части векторов. Наша цель состоит в переформулировке ИУ

(1) и эквивалентного ему ИУ (3) для величин

С и Ф. Поскольку Е - полярный вектор, то

Ф - скаляр, а С - аксиальный вектор (псевдовектор). Выбор этих величин неоднозначен, так как С может быть дополнен градиентом произвольного псевдоскаляра *¥, а потенциал Ф определен с точностью до произвольной константы с0. Для того чтобы исключить неоднозначность в выборе вектора С, подчиним его условию У • С = 0, т.е. будем считать его соленоидальным. Таким образом, введенные величины удовлетворяют соотношениям

у- Е(г ) = -у 2ф ,

Ух Е (г ) = УхУх С =

= ^гай • (Ну -У2)С =-У2С.

(13)

(14)

Далее следует подставить (12) в (3) и произвести разделение вихревых и потенциальных частей на два уравнения. В соотношении (3) сразу выделяем потенциальный вектор

Ео р (г ) = -Уи, где

-о р\

и (г

(г )=- £ С (г - г' )г 1 [Е (г' )У'г](г' +

+ (г -1)£ С (г - г' )к(г' )Е "(г' ) йБ

(15)

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

2 (г -1)£С(г - г')[у'х С(г') - У'Ф(г')](г ' =

V

= к2 (г -1)|£С(г - г')у(г')х С(г')йБ' -

-£УС(г - г')хС(г')(У' - (16)

V

- £у(г')с(г - г')ф(г')(Б'+

+

£Ф(г')УС(г - г')^'1.

При выполнении данного преобразования были использованы теорема о роторе, теорема о градиенте и соотношения векторной алгебры для векторного оператора У' [13]. Возьмем дивергенцию от обеих частей (3). При г £ Б имеем У - Е, поэтому

У • £ С (г - г' )и(г' )х С (г' )(Б' =

Б

= У • £ С (г - г' У(г' )ф(г' )(Б' = -р(ї), г £ Б.

Б

На поверхности указанные величины терпят скачок. В частности, дивергенция второго интеграла равна нулю. Ротор от соотношения (16) в силу (3) совпадает с Ух Е. Поверхностные интегралы разделим на со-леноидальные и потенциальные части:

£ С(г - г' )у(ґ )х С (г' )йБ' = М5 (г) + Мр (г),

Б

£ С (г - г' ) V (г' )ф(г' )(Б' = Р (г)+ Рр (г).

Б

Разделяя соленоидальные и потенциальные части в (3), получаем два уравнения:

Ух С (г ) = к2 (г- 1){М (г )-р (г) +

+ £ С(г )хУ'С(г - г )(V', -УФ(г ) = -к2 (г- 1){Рр (г)-Мр (г) + + у£ Ф (г' ) С (г - г' ) (V' I - У и (г).

(17)

(18)

Общий интеграл последнего уравнения

есть

к

Б

Б

V

ф(Г )= k2 (є-l)x

x |wo (r )-Ф 0 (f)+ J Ф(Г' G (f - r' )dV 'j + (19)

+ U(r)+ Со,

где

Pp(r) = ^Фо(r), Mp(r) = -Wo(r)• (20)

Неизвестные потенциалы Ф0 и W0 удовлетворяют уравнению Пуассона:

V^0 (f) = -V • I v(f' )G(f - f' )ф(Г') dS' =

= I Ф(Г' )V(f' )V' G (f - f') dS' = -p(f),

S

- V • £G (f - f')V(f')x C(f ')dS' =

= -p (f )= V2 W o (f).

(21)

(22)

Решение этого уравнения методом Фурье (методом ФГ) имеет вид

Л

, 1 ^ I dk I dr (2п)3 J J

Фо (f

V^o (f )J rr

p(f' )exp(- jk(f - f'))

(2З)

В (23) первые два интеграла берутся по бесконечной области пространства г' и пространственной спектральной области к = л0кх + + у0ку + г0кг. Используя свернутый вид ФГ

[11], можно (23) представить как

ф„ (?)=

=£ йг'в0 (г - г7)£ ф(г ,)Р(г'Vв(г' - г')йБ*.(24)

Б

Проинтегрируем далее уравнение (17): Vх С(г ) = к2 (е-1) {м, (г)-^ (г) +

+

I C(f ')xV' G(f - f ')dV' [,

C(f) = ko2 (є - l){Cl (f) - C2 (f) + * (f)] . (25) Здесь

Vx C1(f ) = Ms (f) , Vx C2 (f ) = Ps (f) , (26)

С3 (г) = £ С(ї - г')С(г' )(V/ . (27)

V

Для удовлетворения вектором (25) условию У - С = 0 достаточно потребовать, чтобы векторы Сі (і = 1, 2, 3) были соленоидальными.

Поэтому, взяв ротор от соотношений (26) и

(27), получим, что эти неизвестные векторы должны также удовлетворять уравнениям Пуассона с правыми частями - Мя (г),

- Р (г) и - N (г) соответственно. Тогда выражения для этих векторов будут иметь форму (23), в которой для вектора С1 следует поверхностный интеграл из соотношения (21) заменить на ротор вектора М5, для вектора С2 - соответственно на ротор вектора Р, а для С3 - на ротор вектора (27).

В результате получены новые связанные ИУ (19) и (25). Связь осуществляется в силу того, что нормальная компонента поля на Б (поверхностная плотность заряда 7 (г) =

= (г- 1)Е;(г)) в (15) и (19) определяется

потенциальной и соленоидальной частями Е . При к ^ 0 функция (19) определяет потенциал плотности заряда поляризации тела, который согласно принципам электростатики на бесконечности должен стремиться к нулю. Для этого необходимо положить с0 = 0, при этом Ф(г ) = и (г). Полученные

новые ОИУ являются связанными через поверхностные интегралы уравнениями Фред-гольма второго рода с пониженной по сравнению с (1) сингулярностью ядер.

3. Поля в дальней зоне

Приведенные ИУ и ИДУ следует решать в объеме V резонатора, однако представления (1)-(3) справедливы во всем пространстве. В дальней зоне г >> г', поэтому производными в (1) можно пренебречь, а в (2) и (3)

оставить член, пропорциональный |r — r

Перейдем в сферическую систему координат, в которой скалярную ФГ можно представить в следующем виде: G(r - r')« G(r)ejkr cos(r). Здесь как обычно cos(^) = cos(#) cos И+

S

S

2

k

+ sin(#)sin(#')cos(p — <р ) - косинус угла между векторами r иг'. Обозначим En и (0n = (o'n + jo? затухающее во времени собственное поле и комплексную частоту моды с номером n , а>"п > 0. Тогда для собственной

моды имеем

En (? )=°>p(— jfnr / с) F м . (28)

4nrc

Функция F зависит от распределения электрического поля в объеме V и определяет диаграммы направленности излучения ДР для указанной моды [14]:

F {в,р) = j в°cos ^ с [е(г 0— $ (г' )/2 х

V

х sin (e')dr'dO'dp'.

Как видно, поля собственных мод (28) возрастают на бесконечности [15], что связано с экспоненциальным их убыванием во времени: более дальним расстояниям соответствуют более ранние моменты времени излучения (определяемые запаздыванием r/c), когда энергия колебаний в объеме V была экспоненциально больше. Модуль функции

(28) имеет минимум в точке rn = с / (о"п . Если эта точка соответствует дальней зоне, то в заданном направлении внутри V поле может иметь колебательный характер, затем по мере удаления от ДР его модуль убывает до указанного радиуса, а затем начинает возрастать. Для высокодобротных колебаний rn >> a, что оправдывает приближенный метод вычисления действительных собственных частот [1] путем замены полей (28) на убывающие.

4. Численные алгоритмы и результаты

В общем виде все колебания ДР определяются ИДУ (2). Для изотропного диэлектрика можно также использовать ИУ (3), которое для ЦДР имеет вид (5) -(7). Последние уравнения описывают все типы колебаний в изотропном случае. Для однородного ЦДР Ve = 0, поэтому вторые объемные интегралы в этих соотношениях исчезают, а перед первыми получаем спектральный параметр k2 (е — 1). ИУ (5) -(7) можно решить

численно, используя, например, метод обобщенных взвешенных невязок [16] и вытекающие из него методы: Ритца, Трефтца, моментов, Галеркина, коллокаций. В силу комплексности скалярной ФГ спектральный параметр становится комплексным даже в отсутствии потерь в диэлектрике. Общее численное решение спектральной задачи хотя и является строгим и универсальным, но требует нахождения комплексных корней СОп = бо'п + ]&'' характеристического уравнения, обычно представляющего собой определитель с матричными элементами в виде несобственных интегралов в спектральной пространственной области от трансцендентных специальных функций. Необходимость приближенного численного вычисления матричных элементов, усечение бесконечного определителя, ошибки округления при его численном вычислении, приближенное определение его корней может приводить к большим погрешностям, особенно при определении резонансных частот колебаний высших типов. Поэтому целесообразно использовать методы разложения интегрального оператора по малым параметрам. Пусть ДР характеризуется некоторым размером а. В качестве такового могут служить максимальный радиус ЦДР (при этом предположим, что его высота к меньше а). Тогда для низших мод при больших е малым безразмерным спектральным параметром будет к = ка = 2па / X (поскольку Хе~П2 ~ а).

Исходным для построения основанных на этом численных методов служит разложение экспоненты в ФГ (4) в ряд:

в(г - г ') = £(- ]к )пвп (Я ) =

п=0 (29)

= у (- ]кя)п

4пя п! '

Ограничимся N членами такого разложения. В силу равномерной непрерывности (29) в V можно ИУ (3) представить в операторном виде

N N

Е = к2У (- Д)п ЁпЕ -у (- ]к)п £пЕ +

п-° (30)

+у (- Д ) БпЕ,

п=0

где входящие в (30) линейные интегральные операторы выражаются так:

К пЕ = тЧ £ Яп-1 (г - г )[е(г')- 1]Е (г' ,

4т!'1

ЁпЕ = -Ц £ е-1 (г')[Е(г'Ке(г')]х

4тг\ г XV'Яп-1 (г - г')^

Е = -Ц£[е(Г)- 1](Г(г')Е(Г)) 4пп\ -1

(г ))х

XV'Яп-1 (г - г')йБ'

Введем скалярное произведение в пространстве искомых решений:

(и, V )= £ и * (г )у (г )йV .

V

Тогда согласно методу обобщенных взвешенных невязок [16] исходная краевая задача эквивалентна отысканию нулевых стационарных значений (Е = 0) функционала

е = [ Е, Е) - к2у (- )к)п [ Е, т^пЕ) +

п=0

(31)

+у (- jk )п Г Е, £Е) -у (- д )п (Е, ^€,Е').

В силу структуры разложения (29) соотношение (31) является разложением по малому параметру ка. В случае метода Галеркина

Е = Е , и функционал становится квадратичным. Рассмотрим низшие порядки по волновому числу к. При N=0 будем иметь

к 2 = (Е, Е К(Е, £„,£ )-(Е, Б0 Е) (32)

(Еде) '

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

известны. Члены (Е, БгпЕ) с поверхностными и объемными интегралами исчезают только для азимутально-симметричных Н-колеба-

ний, а члены (Е, ЗпЕ) - в случае кусочнопостоянной проницаемости. Наличие поверхностного интеграла свидетельствует о связанной плотности поверхностного заряда, возникающей из-за поляризации диэлектрика. Для высокодобротных колебаний этот интеграл должен либо отсутствовать, либо быть малым, поскольку излучение главным образом связано с соответствующим ему дипольным током поляризации. Разложим поле по полной ортонормированной в V системе векторных функций (в качестве таковых можно взять, например, объемные кусочнопостоянные векторные конечные элементы (КЭ)):

Е(Г ) = £ Етвт (Г ) .

(33)

т=1

В случае объемных КЭ в каждом элементе имеется три вектора с тремя ортогональными поляризациями, поэтому число разбиений объема V равно М/3. Разложение (33) приводит к уравнению в виде равенства нулю определителя размерности М X М , в каждом матричном элементе которого имеется член

7 N+2 гр <—

к максимальной степени. Таким образом, характеристическое уравнение представляет собой полином по к степени М ^+ 2) с действительными и мнимыми коэффициентами, который можно записать как

д(к) = D0 - jK1 - к2+ jК3В3 +

- jMN кМ(N+2)

= ^«1(д,т )= 0.

+ к*Д -jMNк“^N+2>Dм(N+2)- (34)

Здесь введены безразмерные матричные элементы, которые имеют вид

N

Атт' = 8тт + X (- jK+

п=0

N N+2

+у (- jк)nc(:m, = у (- jк)nd(mm,,

п=0 п=0

е(п\ = а_(п+2% К г )

тт т п т

с:: = а -п {(гт , Ёпе.т, )- (ет , )},

при этом йтт' = 8тт'8п0 + Сттт для п < 2,

йн = е тт2)+стт, для 2 < п < N и

V

п

Б

п=0

п=0

= е(тП-Р при п > N. Входящие в (34) коэффициенты выражаются через определители, в которых учтены только матричные элементы, дающие указанную степень спектрального параметра. Например,

д=ае^,2), о№+2)= дм(<С2>).

Нахождение корней полинома является более простой задачей, нежели поиск комплексных корней определителя с матричными элементами в виде интегралов от трансцендентных функций спектрального параметра. Однако порядок многочлена (34) может быть очень большим, поэтому для его уменьшения необходимо согласовывать параметры N и М. Если необходимо найти первые К собственных значений, то целесообразно выбрать К = N = М. Если же для некоторой моды из физических соображений приближенно известна собственная функция Е, то уравнение (31) приобретает вид полнома Б(к ) = Е = 0, где величина Е (31) определена

при Е = Е . Ограничившись N = 2, получим уравнение четвертой степени, решение которого известно аналитически.

Рассмотрим другой подход, основанный на асимптотическом [17] разложении ФГ типа (9) по малому параметру ка. Для колебаний с азимутальным индексом т скалярную ФГ можно записать так:

0(р, р\ф-ф , г - г' ) = _ехр(- ]т(ф-ф))

4п

х

(35)

х {- ;о(1) (р р , г - г') + 0<2) (р р , г - г')}. Здесь первая функция комплексная и имеет

вид

G(l)(р,р/, г - г' ) =

к

=I

1.

Ьр)3т Хррехр- ^к2 -X2 \г - )

2 2 2 X

XX,

(36)

а вторая - действительная и представляется интегралом

0(2)(р, р , г - г') =

_} 1т (2р)3т иррехр)-^2 - к2

г - г

XX

(37)

к л/х2 - к 2

В интеграле (36) сделаем замену X = ка, а в интеграле (37) разложим квадратные корни по малому параметру (к / X)2 < 1:

(X2 - к2 )-1/2 = X 1Ё^# (к / X)’* • п=0 (2п)!!

X2 - к > г=^1 - Ё

(2п - 3)!!

■(к / х)2n^. (38)

Поскольку ряды (38) сходятся неравномерно, соответствующие разложения следует трактовать как асимптотические [17]. Первая функция приобретает вид

0{1\р,р, г - г ) =

1

=к 11т (кар)1т (кар)х

0

{е о^кл/ 1-а2 (г - г'))- ] 81п(кл/ 1-а2 |г - г'!)}

х

лЯ-

айа

а

и при к ^ 0 стремится к нулю, как к2т+1. Ее разложение по параметру к дается разложением в ряды функций Бесселя и тригонометрических функций. Оставляя члены не более второго порядка (что дает четвертый порядок по к), получим:

О

(1)

к 2т+1рр 1

а

к 2а2 т+3

0 I (т

(т!)2-\/1 - а2

л/Т-

гХ

а

х

(р1 + р'2) , ((1 -аа )(г - г')г)

4(т +1)!

2(т!)

Да

2т+1 I

г - г

(т!)2

йа.

Вторая функция выражается в форме

011 ’ = Ё 1т (ар)1т (XP')eХ

п=0 (2п)«! к

Х<

к 2п / к 2п+2 /^2п+4

- +

-+■

2X2 п+1 8х2 п+3

+...

|г - г 1 + ..Л dх.

Входящие в нее интегралы последовательным интегрированием по частям разлагаются

п

2

2

т

0

в ряд по параметру к . Например, после однократного интегрирования

:kn | x~njm {%p)Jm ixpy

~X\Z _ Z

dX:

k

" Jm (kp )jm (kp')e

-k\z_ z'\

n _ 1

^ ЛжІ-n ■

kn I ^—[/» (xp)J. (xpO'

t 1 - n

_X| Z _ z I

\dX-

Однако такая процедура неудобна. Поэтому сделаем в (18) замену а = -\]х2 - к2 и разложим в ряд Тейлора в окрестности

нуля по параметру к2:

G

(2)

= I J. pa2 + k2 )]т (pVa2 + k2 )x

p(- alz - zl) da = g02) + k 2G{2) + k AG2) +....

x exp

Здесь первые два члена имеют вид

G02) = j Jm {Pa)Jm p'a)exp(- аz - z'\)da ,

о

G2(2) = 2 j a'1 Pm (ра )jm (pa)+

2 0

+ P'Jm (pa )J'm (pa)}exp (- a\z - z'|)da.

Рассмотренные ряды являются асимптотическими, и их необходимо обрывать при достижении заданной точности [17].

Очень удобными для численного анализа ДР являются итерационные методы одновременного решения ИУ (ИДУ) и характеристического уравнения. На рис. 1-5 в качестве примера приведены результаты анализа H 01S

и H011 мод однородного и неоднородного

ЦДР. Использованы метод прямой итерации в форме метода последовательных приближений (МПП) и метод минимальных невязок (ММН) с замораживанием значений ФГ от спектрального параметра на предыдущем шаге, примененные к ИУ и строгому характеристическому уравнению. Были использованы одномерные кусочно-постоянные и дифференцируемые (в виде полиномов второго порядка) КЭ, заданные по трем узлам [18]. Оба метода сходятся к одним и тем же результатам за несколько итераций по параметру к.

Объемные (двумерные) КЭ строились в виде прямого произведения одномерных КЭ. На каждой итерации по параметру к вычислялась зависящая от него матрица и выполнялось несколько (или несколько десятков) итераций для получения решения ИУ в виде собственной функции с нормировкой на единицу. Рис. 1 демонстрирует сходимость в за-

f

N

Рис. 1. Сходимость результатов для собственных частот / (ГГц) ЦДР от числа базисных функций N при кусочнопостоянной (1, 2) и полиномиальной квадратичной (3) аппроксимациях для £ = 38, г0 = 5 мм: к = 7 мм (1); к = 4 мм (2, 3)

1, 6/10

£

Рис. 2. Зависимости резонансной частоты/ (ГГц) и добротности 6 однородного ЦДР г0 = 5 мм, к = 7 мм от значения диэлектрической проницаемости: 1 - модаН0\е, 2 - модаН0ц

k

0

0 / к

Рис. 3. Зависимость резонансной частоты / (ГГц) от формы ЦДР для однородного (1, 2) диэлектрика с £= 50 и неоднородного вдоль оси г диэлектрика (3,4) для к = 5 мм: 1 - модаН0\е; 2, 3,4 - модаН011

Е»

Р

Рис. 4. Зависимость Ке(£ф) (сплошные кривые) и 1т(£ф) (штриховые кривые) от координаты р (см) для ЦДР £ = 100, г0 = к = 5 мм: 1 - мода Н01е при г = 0.09; 2 - модаН011 при г = 2.41; 3 - мода Н011 при г = 0.09

висимости от числа разбиений (узлов) по каждой из двух координат. Число пробных функций (порядок матрицы) здесь N . При использовании гладких КЭ размерность задачи снижается примерно в 4 раза. Результаты для собственных частот примерно на 1-1.5% ниже, чем данные из [4]. В частности,

Рис. 5. Зависимость Re(E,p) (сплошные кривые) и 1ш(Еф) (штриховые кривые) от координаты z (мм) для ЦДР £= 100, r0 = к = 5 мм: 1, 2 - модаН01е при р= 3, 4 и p = 1.66 мм; 2 - мода Н011 при р = 4.46 мм

для кривых (2, 3) уточнение по методу Эйт-кена дает f = 5.237524 ГГц, тогда как в [4] f = 5.289 ГГц. С другой стороны, решение самосопряженной задачи с линейным вхождением собственного значения к -2 (т.е. для ФГ при к = 0) рассмотренными методами дает значения, практически совпадающие с [4]. На рис. 2 представлены резонансные частоты и добротности однородного ЦДР как функции диэлектрической проницаемости. Частоты Н011 колебаний при рассмотренной форме ДР примерно в полтора раза, а добротности - на порядок выше, чем для низшего типа Н01§ колебаний, поскольку в первом случае Q связана с магнитно-дипольным излучением, в во втором - с магнитно-квадру-польным. Рис. 3 представляет зависимость резонансной частоты от формы однородного и неоднородного вдоль оси z диэлектрика. Кривая 3 соответствует неоднородному диэлектрику с проницаемостью ^(z) = 25(1 + + cos(^z / к)), а кривая 4 - проницаемости ^(z) = 1 + 49 cos (nz / к). В последнем случае при z = ± к /2 отсутствует скачок £. На рис. 4 и 5 приведены распределения электрического поля в зависимости от координат p и z.

Для обеих мод по обеим координатам обра-

z

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

Заключение

Получены ИУ и ИДУ, описывающие собственные колебания ДР произвольной формы, которые конкретизированы для ЦДР. Предложен метод асимптотического разложения ядра по малому параметру, позволяющий получать характеристическое уравнение в виде полинома по спектральному параметру. Предложены и реализованы для ЦДР эффективные итерационные алгоритмы решения задач на собственные значения. Задачи о возбуждении ДР следует формулировать добавлением к полученным уравнениям неоднородных членов, соответствующих возбуждающему полю. В реальных конструкциях ДР расположен на элементах крепления, поэтому важно учитывать наличие металлических тел. Обобщение приведенных уравнений при наличии бесконечного идеально проводящего экрана получается методом изображений для скалярной ФГ, что, в частности, использовано для моды Н011. Для двух параллельных экранов необходимо уже учитывать бесконечное число изображений. ДР в прямоугольном волноводе и в прямоугольном экране можно проанализировать с помощью метода двумерных и трехмерных изображений [19, 20]. Соответствующие ряды для ФГ хорошо сходятся. Наличие произвольных экранов или металлических тел требует введения члена в виде поверхностного интеграла от искомого распределения поверхностной плотности тока, что может существенно усложнить анализ, поскольку требует решения комбинированного поверхностно-объемного ИУ для распределения поля и тока. Тем не менее рассмотренные методы представляются более общими и перспективными для численных расчетов, чем МЧО и связанный с ним метод ПИУ, поскольку последние применимы только к ДР с координатными поверхностями, а метод ПИУ [3] - к однородным ДР.

Библиографический список

1. Ильченко М.Е., Взятышев В.Ф., Гассанов Л.Г., Безбородов Ю.М., Бергер М.Н., Добромыслов В.С., Капилевич Б.Ю., Нарытник Г.Н., Федоров В.Б., Черний Б.С. Диэлектрические резонаторы. М.: Радио и связь, 1989. 328 с.

2. Van Bladel J. On the resonances of dielectric resonator of very high permittivity // IEEE Trans. Vol.MTT-23, №2. P.199-208.

3. Васильев Е.Н. Возбуждение тел вращения. М.: Радио и связь, 1987. 272 с.

4. Гольдберг Л.Б., Пензяков В.В. Расчет аксиально-симметричных H-колебаний в диэлектрических резонаторах методом интегральных уравнений // РЭ. 1982. Т.27, №9. С.1735-1740.

5. Эскин Г.Э. Краевые задачи для эллиптических псевдо-дифференциальных уравнений. М.: Наука. 1973. 232 с.

6. Тейлор М. Псевдодифференциальные операторы. М.: Мир, 1985. 472 с.

7. Ильинский А.С., Смирнов Ю.Г. Дифракция электромагнитных волн на проводящих экранах (псевдодифференци-альные операторы в задачах дифракции). М.: ИПРЖР, 1996. 176 с.

8. Davidovich M.V. Regularization of kernels for volume integral equations // Modeling in applied electromagnetics and electronics. Saratov: Saratov University Press, 2006. Iss.7. P.30-38.

9. Вычислительные методы в электродинамике / Под ред. Р. Миттры. М.: Мир, 1977. 486 с.

10. Гюнтер Н.М. Теория потенциала и ее применения к основным задачам математической физики. М.: ГИТТЛ, 1953. 416 с.

11. Марков Г.Т., Чаплин А.Ф. Возбуждение электромагнитных волн. М.: Радио и связь, 1983. 296 с.

12. Дао Л.Ц. Математические методы в физике. М.: Мир, 1965. 296 с.

13. Корн Г., Корн Т. Справочник по математике. М.: Наука, 1973. 832 с.

14. Гольдштейн Л.Д., Зернов Н.В. Электромагнитные волны. М.: Сов. радио, 1971. 664 с.

15. Вайнштейн Л.А. Открытые резонаторы и открытые волноводы. М.: Сов. радио, 1966. 476 с.

16. Бреббия К., Теллес Ж., Вроубел Л. Методы граничных элементов. М.: Мир, 1987. 524 с.

17. Олвер Ф. Асимптотика и специальные функции. М.: Наука, 1990. 528 с.

18. Давидович М.В. Метод конечных элементов в пространственно-временной области для нестационарной электродинамики // ЖТФ. 2006. Т.76, вып.1. С.13-23.

19. Давидович М.В. К нестационарной теории возбуждения резонатора // РЭ. 2001. Т.46, №10. С.1198-1205.

20. Давидович М.В. К нестационарной теории возбуждения волноводов // Там же. №11. С.1285-1292.