Применение векторного метода конечных элементов для анализа электромагнитного поля в согласованных пленочных СВЧ-резисторах Текст научной статьи по специальности «Физика»

======================================Известия вузов России. Радиоэлектроника. 2003. Вып. 3

удк 519.68: 538.3+538.5

Ю. Г. Соловейчик, М. Э. Рояк, Е. Б. Корытный, В. П. Разинкин

Новосибирский государственный технический университет

Применение векторного метода конечных элементов для анализа электромагнитного поля в согласованных пленочных СВЧ-резисторах*

Рассматриваются возможности использования векторного метода конечных элементов для численного моделирования трехмерных электромагнитных полей в согласованных пленочных СВЧ-резисторах. Обсуждается вычислительная схема конечно-элементного моделирования с использованием тетраэдральных сеток. Приводятся результаты моделирования одного из вариантов пленочного резистора для частоты 5 ГГц.

СВЧ-Устройства, численное моделирование, электромагнитное поле, метод конечных элементов

В СВЧ-диапазоне для построения мощных широкополосных нагрузок применяют пленочные резисторы сосредоточенного или распределенного типа. В пленочных резисторах мощностью 100 Вт и более на невысоких частотах (порядка сотен мегагерц) ток, протекающий по резистивной пленке, практически не имеет пространственной дифференциации из-за относительно невысокой удельной проводимости пленки и пренебрежимо малых токов смещения. Поэтому при определении входного импеданса и источников тепловыделения для моделирования температурных полей достаточно использовать относительно простые двухмерные математические модели электромагнитного поля. Для частот же порядка 1 ГГц и выше ток в пленке изменяется не только (и даже не столько) в ее поперечном сечении, но и вдоль нее (из-за существенного влияния токов смещения). В этом случае достаточно точные оценки характеристик резистора можно получить лишь на основе трехмерного моделирования электромагнитного поля.

В настоящей работе продемонстрированы возможности использования векторного метода конечных элементов (МКЭ) для численного моделирования трехмерных электромагнитных полей СВЧ-резисторов на частотах, превышающих 1 ГГц.

Постановка задачи. Математическое моделирование электромагнитных процессов будем проводить для конструкции (рис. 1), состоящей из резистивной пленки размером

о о

6 х 8 х 0.006 мм с удельной проводимостью 5 = 4.17 -10 См/М и двух медных контактов

* Статья предоставлена Восточной региональной секцией редакционного совета журнала. © Ю. Г. Соловейчик, М. Э. Рояк, Е. Б. Корытный, В. П. Разинкин, 2003 71

Известия вузов России. Радиоэлектроника. 2003. Вып. 3= z¡

\

X \

Медный контакт

\Резистивная пленка

Л_

Диэлектрическая подложка

Идеальный проводник

Медный контакт

Сторонний ток

\ Идеально проводящее основание \

ч.__________________________________________\

Рис. 1

3 3

размерами 6 х10 х 0.006 мм и 6 х1х 0.006 мм , нанесенных на диэлектрическую подлож-

3

ку размером 10 х19 х 4 мм с диэлектрической проницаемостью 8 = 680 (S0 - диэлектрическая проницаемость вакуума). Сторонний ток, имитирующий генератор, задан в отдельной подобласти размером 1х 0.1 х 4 мм , соединяющей медный контакт с идеально проводящим основанием. Второй медный контакт соединен с идеальным проводником, соединяющим резистор с идеально проводящим основанием. Вся конструкция закрыта идеально проводящим экраном, удаленным от резистора на 60...70 мм. Поскольку задача моделирования электромагнитного поля имеет плоскость симметрии х = 0, разрезающую конструкцию на две одинаковые части, в расчетную область может быть включена только половина описанной конструкции.

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

rotH = JCT +ctE + s(dE/dt); rotE = -8B/dt; divB = 0, (1)

где H - напряженность магнитного поля; J ст - вектор плотностей сторонних токов (возбуждающих электромагнитное поле); а - удельная проводимость среды; E - напряженность электрического поля; t - время; B = цН - индукция магнитного поля (ц - магнитная проницаемость).

Введение вектор-потенциала A, определяемого соотношениями B = rot A; E = -cA/ dt, позволяет получить все требуемые характеристики электромагнитного поля решением только одного векторного уравнения:

rot

ц 1rotA) + а (dA/ dt) + в (d2 a/ dt2) = J ст.

(2)

Обратим внимание на то, что вектор-потенциал А, удовлетворяющий уравнению (2), должен быть разрывным на границах разрыва а и в, т. е. на границах сред с различными

======================================Известия вузов России. Радиоэлектроника. 2003. Вып. 3

удельной проводимостью и диэлектрической проницаемостью. Таким образом, решение этого уравнения может быть найдено только в пространстве разрывных функций, для которых, во-первых, допустимо применение операции rot, а, во-вторых, нормальные составляющие rot A являются непрерывными (поскольку вектор B = rot A должен иметь непрерывные нормальные составляющие вследствие уравнения (1)). Такими свойствами обладают специальные вектор-функции, используемые в векторном МКЭ [1]. Поэтому в данной статье рассматривается вычислительная схема решения поставленной задачи с применением указанного метода.

Вариационная постановка. Аппроксимация по времени (2) приводит к векторному уравнению

rot (ц- 1rotA ) + уА = F. (3 )

В случае, когда сторонние токи изменяются во времени по гармоническому закону с круговой частотой ю, входящие в уравнение (3) величины A, F и у являются комплекс-

2

ными (F = JCT, у = jam - вю , где j - мнимая единица). В общем случае, когда сторонние токи являются произвольной функцией времени, все входящие в уравнение (3) величины вещественны, причем коэффициент у и вектор-функция F определяются разностной схемой аппроксимации по времени. Например, для трехслойной полностью неявной схемы с постоянным шагом At коэффициент у и вектор-функция F имеют вид:

3^ £ ^ т

y=—+—т; f = jct-

2a (f ст

a

2Д t

( At )2

A-2 +

2a 2s

+ -

& (At)2

A -1,

A —2 4-1

где A и A - значения вектор-потенциала на предыдущих временных слоях, а решение A уравнения (3) является значением вектор-потенциала на текущем временном слое.

Получим вариационную формулировку для уравнения (3) в форме Галеркина. Для этого домножим векторное уравнение (3) скалярно на пробную вектор-функцию ¥ и проинтегрируем полученное скалярное уравнение по расчетной области Q :

J rot (^-1rotA ) WQ + J yAWQ = J F WQ . (4)

Q Q Q

Применим к первому слагаемому уравнения (4) векторную формулу Грина:

J ^-1rotArotWQ - JV-1 [( rotA) x n ] WS + J yAWQ = J F WQ , (5)

Q S Q Q

где S - граница Q; n - внешняя по отношению к Q нормаль к S, а символ " х" обозначает векторное произведение.

На тех частях Sp границы S, которые соответствуют поверхности идеального проводника, должно быть задано краевое условие E х n = 0, т. е. касательные составляющие E на поверхности идеального проводника должны равняться нулю. Соответственно, на

Известия вузов России. Радиоэлектроника. 2003. Вып. 3======================================

Se должно выполняться аналогичное уравнение и для вектор-потенциала A и пробных

вектор-функций ¥:

A х n = 0, Y х n = 0. (6)

С учетом свойств смешанного произведения векторов и соотношения (6) получим [( rot A) x n ] Y = (rot A)(n x Y) = ( rot A) 0 = 0,

т. е. на границах Sp поверхностный интеграл в уравнении (5) равен нулю.

Рассмотрим поверхность Sв, на которой задано условие симметрии B х n = 0 (т. е.

на Sb ненулевыми являются только нормальные составляющие поля B). Тогда и на Sв

[(rot A) x n] Y = (B x n) Y = 0Y = 0.

Таким образом, в рассматриваемой задаче подынтегральное выражение в поверхностном интеграле уравнения (5) равно нулю на всей границе S расчетной области Q и вариационная постановка принимает окончательный вид:

J ц-1 rot A rot WQ + J yAWQ = J F WQ. (7)

Q Q Q

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

p (i, j) соответствия пар глобальных номеров вершин сетки i и j глобальным номерам ее ребер. На каждом из тетраэдров Q^ , содержащих узлы сетки i и j (и, соответственно, ребро p(i, j)), базисная вектор-функция Yj) будет иметь вид [1]:

ТР(Ü) = Li grad Lj + Lj grad Li, где i и j - локальные номера вершин i и j, а Li и Lj - L-координаты рассматриваемого тетраэдра Q^ (называемые также барицентрическими, или естественными координатами тетраэдра) [2]. Под L-координатами на тетраэдре Q^ понимаются линейные функции координат вида

L¡ ( X, y, z ) = а0 + a^x + а2y + a3 z, (8)

такие, что Li равна единице в вершине тетраэдра с номером i и нулю - во всех остальных

его вершинах. Коэффициенты a'¡ в выражении (8) могут быть найдены из соотношения

Г 1 1 1 О

ССд СС1 СС2 СС3

2 2 2 2 ССд СС1 СС2 СС3 _1

3 3 3 3 ССд СС1 СС2 СС3

4 4 4 4 ^а о а1 а 2 щ

О

(9)

где матрица Ъ определяется через координаты (х, у, ) 1-й вершины тетраэдра следующим образом:

(10)

'1 1 1 о

Х1 Х2 Х3 Х4

У1 У2 У3 У 4 ^ ^2 ^3 ^4)

Заметим, что вектор-функции ¥^ л) имеют непрерывные касательные составляющие на всех гранях тетраэдров, а разрывы их нормальных составляющих позволяют естественным образом учитывать разрывы вектор-потенциала А на поверхностях, где могут быть разрывны коэффициенты а и 8 .

Будем искать решение А уравнения (7) в виде разложения по базисным функциям ¥т,

р (г,]): А = ^ , где цт - искомые коэффициенты разложения. Заменяя пооче-

т =

редно пробную функцию ¥ на базисные вектор-функции, получим систему уравнений

^ Чт |^-1Г01 ТтГС1 + ^ Чт |У *т*I|.

т О т о о

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

Оя = Г, (11)

в которой элементы матрицы О и вектора Г определяются соотношениями:

От/ = |ц-1го1 ТтГ01+ |у ¥т; Г/ = \Р^О .

О О □

Глобальная матрица О СЛАУ (11) может быть вычислена как сумма вкладов (локальных матриц) от каждого тетраэдра О.^ . Формулы для вычисления локальных матриц конечных элементов приводятся, например в работах [3], [4]. Авторами статьи будут использоваться более удобные для программирования формулы:

| ц-1гог Тт гог Т^О = (4/ц)| Qk | (аа х аЬ ас х ай ) ;

| УЧт0 = У(^achd -^adPbc +ЧёPac)

т

Известия вузов России. Радиоэлектроника. 2003. Вып. 3======================================

где m = p (а, b); l = p (c, d) (a, b, c, d - глобальные номера вершин тетраэдра, соответст-

вующие локальным номерам a, b, c, d ); векторы коэффициентов аг = , а2, a^j определяются из соотношения (9), а коэффициенты Ху и Pj - соотношениями [2] :

i i г I 0.1mesQk, i = j;

P;J = grad Z, grad Z ,■ = a aJ ; = Z,Z ,dQ. = < lJ 1 J l] J lJ |0.05mesQk, /*/,

Qk 1 k'

где mes Qk = (1/6) |det D, а матрица D определена соотношением (10), т. е. Ху и Pj -

компоненты локальных матриц массы и жесткости стандартных узловых линейных элементов на тетраэдрах.

Результаты численного моделирования. Численные расчеты проводились по изложенной ранее вычислительной схеме для гармонических сторонних токов с частотами в диапазоне от 200 МГц до 5 ГГц. Результаты, полученные для частот до 500 МГц, практически не отличались от результатов двухмерного моделирования, когда считалось, что токи текут только параллельно продольной оси резистора, а магнитное поле имеет только перпендикулярные к этой оси составляющие. На более же высоких частотах, когда с проводников стекают существенные токи смещения, картина электромагнитного поля может кардинально отличаться от полученной с помощью двухмерного моделирования.

Приведем результаты численного моделирования электромагнитного поля для частоты 5 ГГц. Сторонний ток был задан таким, чтобы выделяемая на резистивной пленке тепловая мощность составила 200 Вт, что соответствует согласованному режиму.

На рис. 2 представлено распределение вектора плотности тока в продольном сечении резистивной пленки, перпендикулярном плоскости симметрии, в моменты времени, когда сторонний ток максимален (а) и равен нулю (б). Из рисунка видно, что поле токов очень дифференцированно не только поперек, но и (и даже более существенно) вдоль пленки.

Это объясняется существенным влиянием токов смещения, распределение вектора плотности которых в плоскости симметрии х = 0 приведено на рис. 3 (рис. 3, а соответствует максимальному стороннему току, а рис. 3, б - току, равному нулю). Распределение вектора поверхностной плотности тока Jn0B в те же моменты времени на идеально проводящем основании представлено на рис. 4. Поверхностная плотность тока на идеальном проводнике определяется соотношением Jпов = H х n, где n - нормаль к поверхности

идеального проводника.

Рис. 5 иллюстрирует распределение плотности объемных источников тепловыделения p ( х, y ), порожденных токами в резистивной пленке. Как видно из этого рисунка, их

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

Плотность тока: -150 мА/мм2 ^ 300 мА/мм2

Рис. 2

z, мм

У У М И У МММ?

И I | | I I

М У У ? ! У У У ?

П^4 1 п , v ' '

q ^ s \ у < •

1 Ч * ' '

У у у í Í ? Í 1 I ' '

. I I I » *

^ I

У У У У У У ммм

ч >1 ^ -«» —

1П|И111ПП11ИЧ"|'Г тгсгт

0

- 11- 10- 9- 8- 7- 6- 5- 4- 3- 2- 1 0 1 2 3 4 5 6 7 8 у, мм

а

z, мм

0

i i é i 4

i i I i I

у V y y

y y

y y

? У 44-

ito m

У У ?

f У »

? У »

У У »

У У '

f У 1У '

11 - 10- 9- 8- 7- 6- 5- 4- 3- 2- 1 0 1 2 3 4 5 6 7 8 У, мм

б

Плотность тока смещения: 100 мкЛ/мм2

Рис. 3

200 мкА/мм

2

4

2

4

2

x, мм ............................

0 - 12

I г I Г I I

1 Ч \ ^

ч

<4

/ /■>>>>>• ^ -/ ^ ^ >> ^ ^ -р >>->»-

- 10 - 8 - 6 - 4 - 2

10 у, мм

а

x, мм

- 12 - 10 - 8 - 6 - 4 - 2

Плотность тока смещения:

0

б

2

10 у, мм

0.75 мкА/ мм Рис. 4

2

1. 5 мкА/ мм

2

p = 2.8 кВт/мм 2.4 1.2 1.0

3

x, мм

0.8 0.4 0.6

2.6 2.0 1.6

1.4

Медный контакт

1.0

2.0

0.4

0.2

0.2

Резистивная пленка

2.4 1.8

1.4

Медный контакт

0

4

Рис. 5

у, мм

5

0

2

4

6

8

5

0

4

6

8

6

2

======================================Известия вузов России. Радиоэлектроника. 2003. Вып. 3

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

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

1. Bossavit A. Whitney forms: a class of finite elements for three-dimensional computations in electromagnet-ism // IEE Proc. 1988. Vol. 135. Pt. A. P. 493-500.

2. Норри Д., Фриз де Ж. Введение в метод конечных элементов. М.: Мир, 1981. 304 с.

3. Lee J-F, Mittra R. A note on the application of edge-elements for modeling three-dimensional inhomogene-ously-filled cavities // IEEE Trans. Microwave Theory and Techniques. 1992. Vol. MTT-40. P. 1767-1773.

4. Tharf M. S., Costache G. I. Finite element method solutions of field distributions in large cavities // International Journal of Numerical Modelling: Electronic Networks, Devices and Fields. 1994. Vol. 7. P. 343-355.

Y. G. Soloveitchik, M. E. Royak, E. B. Korytniy, V. P. Razinkin

Novosibirsk state technical university

The Vector Finite Elements Method Using to Analyze Electromagnetic Field in Coherent SHF Film Resistor

The capability of vector finite elements method using to numerical modeling of three-dimensional electromagnetic fields in coherent SHF film resistor are considered. The computing circuit of finite element modeling using tetrahedron meshes is discussed. The results of modeling of one film resistor variant for 5 GHz frequency are given.

SHF-units; numerical modeling; electromagnetic field; finite element method

Статья поступила в редакцию 7 июля 2003 г.