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

радюф1зика

радиофизика _карюрну8!с8_

УДК 621.372.2:004.942

Л. М. Карпуков, Р. Д. Пулов

ПРОЦЕДУРА КВАЗИСТАТИЧЕСКОГО МОДЕЛИРОВАНИЯ МИКРОПОЛОСКОВЫХ СТРУКТУР

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

ВВЕДЕНИЕ

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

© Карпуков Л. М., Пулов Р. Д., 2005

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

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

Следует отметить, что в практике квазистатических расчетов микрополосковых устройств получили распространение две разновидности интегральных уравнений [3, 4]. Интегральное уравнение, используемое в [3], легко составляется с помощью известной функции Грина для свободного пространства, однако оно имеет широкую область определения, включающую в себя поверхности проводников, экранов, границ раздела диэлектрических сред. Область определения интегрального уравнения, используемого в [4], ограничена только поверхностью проводников. Необходимость минимизации размерности вычислительных задач обусловила выбор интегрального уравнения с компактной обла-

стью определения для рассматриваемом ниже процедуры квазистатического моделирования. Составление этого уравнения требует предварительного определения функции Грина подложки моделируемого микро-полоскового устройства. В случае многослойной подложки нахождение функции Грина представляет собой сложную задачу. Ее решение существенно упрощается при проведении расчетов в спектральной области [5]. Однако переход от спектральных представлений функций Грина к оригиналам с помощью интегрального преобразования Фурье связан с большим объемом вычислений. В [6] предложен альтернативный метод составления явных зависимостей для функций Грина, основанный на разложении их спектральных представлений в ряды по функциям, оригиналы которых известны. Ниже на основе этого метода составлены формулы, обеспечивающие эффективный расчет явных зависимостей для функций Грина многослойных структур. Эти формулы использованы при алгебраизации интегрального уравнения по методу моментов. С их помощью получены простые и экономичные соотношения для расчета коэффициентов матриц систем алгебраических уравнений. Для повышения эффективности решения систем уравнений большой размерности применено вейвлет-преобразование.

МЕТОД РАСЧЕТА ФУНКЦИИ ГРИНА МНОГОСЛОЙНЫХ СТРУКТУР

Для квазистатического моделирования в работе используется интегральное уравнение, связывающее потенциал ф и поверхностную плотность ст заряда на токонесущих проводниках исследуемой конструкции микрополоскового устройства [4]:

ф(г) = |G(г, г')ст(г')dr',

(1)

где G(г, г') - функция Грина структуры подложки микрополоскового устройства, 5 - поверхность проводников.

Эффективность моделирования на основе уравнения (1) во многом определяется методом вычисления функций Грина структур подложек, которые в общем случае могут быть многослойными.

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

Ега+1

Т"?

о 4

<> 3

а)

б)

Рисунок 1

ми проницаемостями 81, 8п + 1 или между электрическими или магнитными экранами, г-й слой структуры характеризуется относительной диэлектрической проницаемостью 8, и толщиной к.

Введем декомпозиционную модель исследуемой структуры и представим ее декомпозиционной схемой, изображенной на рис. 1, б [6].

Составим из матриц рассеяния, моделирующих границы раздела диэлектрических слоев, сводную блочно-диагональную матрицу и запишем в спектральной области систему уравнений для комплексных амплитуд отраженных и- и падающих и + волн с учетом входных воздействий ф, в граничных сечениях, отмеченных на рис. 1, б цифрами [6]. Исключим из уравнений амплитуды падающих волн с помощью условий связи для волн в г-м слое:

»+ и = и

-к,к

г + 1е

ч +1

(2)

В результате получим математическую модель исследуемой структуры:

и = 5 • I • е

-к,к-

и + ТФ,

и = I • е

-М.

и

(3)

(4)

Здесь векторы и , и , Ф составлены соответственно из амплитуд отраженных и падающих волн и входных воздействий в граничных сечениях. Матрицы 5, Т, I имеют блочно-диагональную форму:

5 = diag

-Г

12'

Г

23

1 - Г

23

1 + Г23 -Г23

Г

п - 1,

1 - Г

1 + Г

-1,п

п - 1, п

' п - 1, п

,Г

п, п - 1

, Ггк =

8г + 8 к

1607-3274 «Радюелектрошка. 1нформатика. Управлшня» № 2, 2005

кгк

+

= и е

8

Т = diag

( 1 г \

Тр, Т 23 Т23 Т п -1, п Т п -1, п Т ' п, п + 1

V } 23 Т 23 Т п -1, п Т п -1, п У

Тгк =

0 1 1 0

0 1 1 0

0 1 1 0

0 1 1 0

I = diag^

Для трехмерных краевых задач кг = ^к2. + кУ,

-г + ^

Фурье-изображение потенциала р г е Фф от точечного источника в точке воздействия записывается в виде

рр1 = ехр()к.х + }ку)/(4п ), (5)

где 80 - диэлектрическая проницаемость свободного пространства, / - мнимая единица.

Для двумерных задач, например, на плоскости у0г:

рРг = ехР(/кх)/( 2 П802кг), кг = ку.

Представив решение уравнения (3) в виде ряда по степеням ехр(-кгк) и объединив его с выражением (4), составим формулу для расчета потенциалов, определяющих спектральные функции Грина в граничных сечениях:

" ® п -пк,к" п -(п +1)кк"

и = £(5 • I)пТ • е г Ф + I •(Б • I) Т • е Ф.

п = 0

(7)

Осуществив переход на основании данной формулы от спектральных представлений в область оригиналов, получим явную формулу для расчета функций Грина в слоистой среде:

ад

и = £(5 • I)пТ -Ф(п) +1 •(Б • I)пТ -Ф(п + 1), (8)

п=0

где рг е Ф записываются, в соответствии с (5), (6), как

рг (п) = -- 1 , г = л/а X + А у2 (9)

4тс80л/ г + (пк)

рг( п) = г-1-1пл/ Ау2 + (пк)2

2 7С80

(10)

Здесь Ах, Ау - расстояния между точкой возбуждения и точкой наблюдения по соответствующим осям.

Составление формул вида (8) по декомпозиционным схемам анализируемых структур не вызывает за-

труднений. Например, для однослойной структуры при размещении металлического экрана в сечении г\ и при источнике поля в сечении ¿2 (см. рис. 1) функции Грина определяются формулой:

=Т

23 X

п=0

0 -1 Г23 0

0

Р2( п)

0 1 0 -1 п 0

1 0 _Г23 _р2( п + 1)

(11)

При добавлении еще одного слоя:

Т23 х ^

п=0

0 0

р2( п) + I • Бяг • р2( п + 1)

0 0

_ 0 _ _ 0 _

(12)

(6) где 1: =

0 10 0 10 0 0 0 0 0 1 0 0 10

, S1 = Б • I =

23

1 + Г23 0 0

-10 0 0 0 1 - Г23

г

23

0

0 Г34 0

В аналогичной форме записываются соотношения для произвольного числа слоев. Для сокращения вычислительных затрат на расчет по этим соотношениям введем для матрицы Sl преобразование подобия [7]. В результате получим следующую расчетную формулу:

ад

и = Тд X [Е • рч( п) +1 •рч(п + 1)]- М •Хп • Вя. (13)

п=0

Здесь Е - единичная матрица; Х - диагональная матрица, составленная из собственных чисел матрицы М - модальная матрица, составленная из собственных векторов матрицы Sl; Тд - коэффициент Тгк, определенный для сечения д, где приложено возбуждение р9; Вд - столбец матрицы М-1, соответствующий д-му сечению.

МЕТОДИКА РЕШЕНИЯ ИНТЕГРАЛЬНОГО УРАВНЕНИЯ

Алгебраизацию интегрального уравнения (1) произведем методом моментов с представлением результата в виде матричного уравнения [4]:

О^ а = и,

(14)

где

Ор, д = | шр(г)| О(г> г')г')dr,

Р

и

+

и

и

и

2

и

3

и

4

или

Up = J wp( r )ф( r)dr•

Для сокращения объема вычислений в качестве базисных Юц и весовых Шр функций применим кусочно-постоянные функции при равномерном разбиении поверхности проводников. Расчет коэффициентов Gpц уравнения (14), соответствующих функции Грина в сечении г = гр структуры подложки на рис. 1, выполним по соотношению (13) с подстановкой вместо фц значений интегралов

Gq(n) = jj9q(n, r, r')drdr',

(15)

полученных в аналитическом виде от подынтегральных функций (9), (10).

Матрица уравнения (14) является полностью заполненной, поэтому вычислительная стоимость решения этого уравнения высока, она пропорциональна 0(Ы3), N - число уравнений.

Эффективным средством снижения вычислительных затрат на решение больших систем уравнений является вейвлет-преобразование, обеспечивающее получение разреженных матриц, стоимость операций с которыми не превышает 0^2) [8-10].

Уравнение (14) с помощью вейвлет-преобразования записывается в виде:

G

W

tw uW,

(16)

где GW = W • G • , = W •ст, иш = W • и, W -матрица вейвлет-преобразования, Т - знак транспонирования. Решение уравнения определяется как

а = W

т

W

(17)

Для построения матрицы W были использованы вейвлеты Добеши [10]. Генерация вейвлетов осуществлялась с помощью компьютерных программ из [11]. Для примера на рис. 2 представлены значения коэффициентов матрицы вейвлет-преобразования для вейвлетов Добеши с 512 точками дискретизации и 8 моментами.

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

ЧИСЛЕННЫЕ РЕЗУЛЬТАТЫ

Для проверки эффективности разработанной вычислительной процедуры квазистатического моделирования выполнен расчет параметров микрополосковой линии на однослойной подложке с относительной диэлектрической проницаемостью 8Г = 4 и толщиной к = 1 мм.

400

200

0 200 400

Рисунок 2

В таблице 1 приведены данные о продолжительности расчета эффективной диэлектрической проницаемости еэф линии двумя методами А и В. Метод А изложен в [3] и использован в популярной среди разработчиков СВЧ устройств системе проектирования Microwave Office [12]. Метод В предложен в настоящей работе. Расчеты выполнялись на одном и том же компьютере. Сравнивались вычислительные затраты на составление матриц коэффициентов уравнений и их решение методом LU-разложения. Расчеты выполнены при разбиении полоски микрополосковой линии, шириной w, на 16 частей (Nw = 16). Отметим, что число разбиений Nw определяет число уравнений в методе В (N = Nw). В методе А система уравнений имеет большую размерность из-за необходимости введения дополнительных уравнений относительно вторичных источников на границах раздела диэлектриков. В расчете по методу А половина ширины d подложки была взята равной 5h с числом разбиений Nd = 16, что соответствовало числу уравнений N= 48.

Таблица 1 - Сравнение вычислительной эффективности методов

Методы

А В

сек. 0,056 0,024

w / h = 2 еэф 3,0751 3,0768

w / h = 1 еэф 2,9208 2,9162

Эффективность А/В

2,33

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

В таблице 2 приведены данные о погрешности расчета волнового сопротивлении с и эффективной диэлектрической проницаемости 8эф рассматриваемой мик-рополосковой линии по формулам (16), (17) с помощью вейвлет-преобразования при отбрасывании коэффициентов, величины которых не превышают 0,1% от макси-

10

ISSN 1607-3274 «Радюелектрошка. 1нформатика. Управлшня» № 2, 2005

p

мального коэффициента матриц. За точное решение взяты данные, полученные непосредственно из решения уравнения (14) с помощью метода Еи-разложения.

На рис. 3 изображена зависимость времени решения системы уравнений в секундах от ее размерности N при решении методом вейвлет-преобразования (кривая 1) и методом Еи-разложения (кривая 2).

10000

Рисунок 3

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

ЗАКЛЮЧЕНИЕ

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

риц и снижение вычислительных затрат на решение систем уравнений большой размерности.

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

ПЕРЕЧЕНЬ ССЫЛОК

1. Проектирование интегральных устройств СВЧ: Справочник / Ю. Г. Ефремов, В. В. Конин, Б. Д. Солганик и др. - К.: Тэхника, 1990. - 159 с.

2. Миролюбов Н, Н, и др, Методы расчета электростатических полей. - М.: Высшая школа, 1963. - 415 с.

3. Bazdar M, B,, Djordievic A, R,, Harrington R, F, Evaluation of quasi-static matrix parameters for multiconductor transmission lines using Galerkin's method dielectric // IEEE Trans. Microwave Theory Tech. - 1994. - Vol. 42. -No. 7. - P. 1223-1228.

4. Сильвестер П,, Феррари P, Метод конечных элементов для радиоинженеров и инженеров-электриков: Пер. с англ. - М.: Мир, 1986. - 229 с.

5. Bhat B, Koul K, Unified approach to solve a glass of strip and microstrip transmission lines // IEEE Trans. Microwave Theory Tech. - 1982. - Vol. 30. - No. 5. - P. 679-683.

6. Карпуков Л, М, Построение и анализ декомпозиционных моделей микрополосковых структур // Радиоэлектроника. - 1984. - Т. 27. - № 9. - С. 32-36. (Изв. высш. учеб. заведений).

7. Сигорский В, П, Математический аппарат инженера. -К.: Техшка, 1975. - 768 с.

8. Beylkin G,, Coifman R,, Rokhlin V, Fast Wavelet Transform and Numerical Algorithms. - Commun. Pure Appl. Math. - 1991. - P. 141-183.

9. Wagner R, L,, Chew W, C, A Study of Wavelets for the Solution of Electromagnetic Integral Equations. // IEEE Trans. AP. - 1995. - V. 43. - No. 8 - P. 802-810.

10. Daubechues I, Ten Lectures on Wavelets (CBMS-NSF series in Applied Maths #61). Philadelphia: SIAM, 1992.

11. Press W, H,, Teukolsky S, A,, Vetterling W, T,, Flanne-ry B, P, Numerical Recipes in FORTRAN: The Art of Scientific Computing, Second Edition. - New York: Cambridge Univ. Press, 1992.

12. Разевиг В, А,, Потапов Ю, В,, Курушин А, А, Проектирование СВЧ устройств с помощью Microwave Office. -М.: СОЛОН-Пресс, 2003. - 495 с.

Надшшла 21.02.05 Шсля доробки 8.10.05

Розроблено обчислювальж основи процедури складання i рШення ттегральних р1внянь при квазистатичному моделюванн мiкросмужковых пристрою, реалiзованих на багатошаровш тдкладщ. Запропоновано явт формули для розрахунку функцш Грiна багатошарових тдкладок i коефiцieнтiв матриць систем рiвнянь, що виходять при алгебраiзацi'i ттегральних рiвнянь методом моментiв. Для рШення систем алгебраЧчних рiвнянь застосовано вейвлет-перетворювання. Представлено результати чи-сельних розрахунтв, що тдтверджують ефективтсть розробленоЧ процедури моделювання.

The computational foundations of formation and solving of the integral equations at quasi-static modeling of microstrip devices on a multi-layer substrate are developed. The explicit formulas for computation of Green functions of multi-layer substrates and coefficients of matrices obtained during algebraization of integral equations by a moment method are proposed. The wavelet transformation is utilized for solving systems of the algebraic equations. The results of the numerical calculations proving the efficiency of the developed modeling procedure are submitted.

Таблица 2 - Сравнение точности вычислений

Размерность Погрешность вычисления еэф, % Погрешность вычисления р, %

32 0,0022 0,00001045

64 0,0095 0,00714

128 0,015 0,0029

256 0,022 0,025

512 0,03 0,032

1024 0,075 0,19