Оценка точности вычисления функции Грина в слоистой среде Текст научной статьи по специальности «Математика»

Вычислительные технологии

Том 11, № 5, 2006

ТОЧНОСТИ ВЫЧИСЛЕНИЯ ГРИНА В СЛОИСТОЙ СРЕДЕ

И. А. Конников Санкт-Петербург, Россия e-mail: [email protected]

A technique for computing of the Green function for the Laplace equation in a layered medium is proposed. The technique enables to check the accuracy of computation of the potential around a unitary point charge. Method of the solution of the problem relies on properties of the Bessel and Struve functions. An example of calculation is given. The results are valid for a broad class of problems solved using the Green functions technique.

Введение

Вычисление функции Грина для уравнения Лапласа в слоистой среде — задача, простое и эффективное решение которой требуется для многих технических приложений. Так, например, для расчета емкостей коммутационных проводников микросхем и печатных плат требуется найти функцию Грина в квазистационарном приближении, решив уравнение Лапласа для потенциала электрического поля единичного заряда, а для реализации математического и программного обеспечения систем автоматизированного проектирования (САПР) необходимы библиотека математических моделей конструкции проектируемых устройств и адекватные формы и методы представления моделей.

1. Получение модели

Для многих технических приложений в качестве физической модели среды, где имеется точечный единичный заряд, целесообразно принять слоистую среду, состоящую из произвольного числа плоскопараллельных слоев, расположенных перпендикулярно к оси аппликат. При расчете емкостей слои имеют разные диэлектрические проницаемости и считаются гомогенными, изотропными, непроводящими и неограниченными в азимутальном направлении. При наличии металлического корпуса первый и/или последний слой может быть ограничен идеально проводящей плоскостью. Физические процессы в модели считаются квазистационарными. Потенциал электрического поля ip во всех точках такой модели (кроме точки, где расположен единичный точечный заряд) удовлетворяет уравнению Лапласа V2p = 0. Методика определения функции Грина для уравнения Лапласа в такой модели предложена в [1] и развита в целом ряде работ, в том числе [2, 3]. Решение уравнения Лапласа для потенциала единичного точечного заряда в ¿-м слое такой модели (т. е.

© Институт вычислительных технологий Сибирского отделения Российской академии наук, 2006.

ОЦЕНКА ФУНКЦИИ

функция Грина в слоистой среде), как известно [2-4], дается формулой

О = I /о(Лг)Ф4(Л, г - ¿о(1) 4п£о ]

о

где е0 — абсолютная диэлектрическая проницаемость свободного пространства; 30 — функция Бесселя первого рода нулевого порядка; г = \/(ж — х0)2 + (у — у0)2; х0, у0, г0 — абсцисса, ордината и аппликата точки, где вычисляется поле; х, у, г — абсцисса, ордината и аппликата источника заряда; функция Ф1(Л, г — г0) определяется из граничных условий. Метод получения аналитического выражения для функции Ф1(Л,г — г0) в строгом классическом варианте изложен в [2, 4].

Для Ж-слойной физической модели существует N функций Ф1(Л, г — г0), I = 1, N. Для любого слоя коммутации в плоскости г = г0 имеется своя функция ФДЛ,г — г0) = Ф1 (Л). При К-слойной коммутации необходимо использовать К функций Ф1(Л) для всех слоев коммутации (1 = 1,К). При необходимости учитывать межслойные связи надо использовать еще и К!/[2(К — 2)!] функций Ф^(Л, гу — гг), где — — расстояние между слоями коммутации. Множество функций Ф1(Л) и Ф^ (Л, — гг) для всех слоев коммутации является тем формальным математическим объектом, который содержит всю информацию о конструкции проектируемого устройства, необходимую и достаточную для расчета емкостей проводников. Этот объект целесообразно использовать в качестве математической модели конструкции, например, в САПР микросхем и печатных плат на стадии расчета емкостей коммутационных проводников. Параметрами такой математической модели будут диэлектрические проницаемости и толщины слоев физической модели. С точки зрения программной реализации совокупность таких моделей может быть реализована в виде библиотеки математических моделей.

2. Рабочая форма представления модели

Несобственный интеграл (1) выражается через первообразные только в простейших случаях. Предлагаемые различными авторами способы его приближенного вычисления основаны на использовании формулы Вебера — Липшица

те

[ ^(Лг) ехр(—Лт)ЛЛ = *--, т > 0. (2)

3 уг2 + т2

0

Точность вычисления функции Грина при этом либо оценивается эвристически, либо не рассматривается вообще: оценивается лишь точность представления функции Ф(Л). Редким исключением является работа [4], где при решении основной задачи электроразведки функцию Ф(Л) удается представить рядом Маклорена с последующим его интегрированием; в результате функция Грина получается представленной медленно сходящимся рядом. При некоторых сочетаниях значений проводимостей слоев ряд является знакопеременным и тогда нетрудно оценить погрешность, которая получается при замене ряда его отрезком. В [3] предлагается аппроксимировать функцию Ф(Л) полиномами Лежандра по степеням и = ехр(—Лт), где т — нормирующий множитель, причем трудоемкость вычисления коэффициентов аппроксимации сравнима с трудоемкостью вычисления самого интеграла (1). В [5] предлагается аппроксимация интерполяционным многочленом (линейной комбинацией

функций Чебышева первого рода) с оптимальным выбором узлов на интервале [0,1] также по степеням u = exp(- At). Учитывая простоту процедуры получения коэффициентов интерполяции [6] и известные возможности аппроксимации в базисе функций Чебышева первого рода, предлагаемый в [5] вариант на первый взгляд можно считать предпочтительным. Однако следует учесть следующее.

Точность аппроксимации степенным полиномом по степеням переменной u £ [0,1] на границах интервала аппроксимации имеет особое значение. Погрешность приближения Ф(А) на левой границе указанного интервала проявляется как неточность определения свободного члена аппроксимирующего полинома. Слагаемые полинома со старшими степенями переменной u являются лишь поправкой, учитывающей неоднородность среды; абсолютная величина их суммы не превосходит абсолютной величины свободного члена, и свободный член должен рассчитываться с максимально доступной точностью. Кроме того, в зависимости от вида физической модели функция Ф(А) на одной из границ интервала может сколь угодно близко приближаться к нулю или даже обращаться в нуль, поэтому даже незначительная абсолютная погрешность аппроксимации функции Ф(А) вблизи любой из границ интервала может проявляться как большая относительная погрешность вычисления подынтегральной функции интеграла (1). Влияние этой ошибки на точность вычисления функции Грина может оказаться определяющим, необходимо избежать появления участка с большой погрешностью вычисления подынтегральной функции интеграла (1). Но, как известно, аппроксимация в базисе функций Чебышева первого рода позволяет минимизировать максимальную абсолютную ошибку и абсолютная погрешность аппроксимации (если она отлична от нуля) именно на границах интервала аппроксимации максимальна. Поэтому целесообразно формально расширить указанный интервал и интерполировать на интервале [u[,u]] (при этом u[ < 0, U] > 1) так, чтобы при u = 0 и u =1 погрешность аппроксимации отсутствовала; для u < 0 и u > 1 интерполирующий многочлен Л^) не используется. Тогда, принимая формулы из [6, §20] в качестве исходных, получим

Ф(А) « ЛЫ = Ao/2 + V AuTv( 2u - u[ - u]) , (3)

V u] - u[ /

где Y — старшая степень интерполирующего многочлена;

Y

Au = V Ф(А;)cos (2j + 1)Un, и = 0,Y;

v y + 1 ^ K 7 2Y + 2 j=o

Tv — функция Чебышева первого рода степени и;

In Uj u] + u[ , u] - u[ (2j + 1)n .

А7 =--, u7 = —-LH—--L cos——-, j = 0, Y;

7 T ' 7 2 2 2Y + 2

1 — cos a0 1 + cos a0 n (2Y + 1)n

u[ = -, u] =-, ao = ox , aY = ox .

cos aY — cos a0 cos a0 — cos aY 2Y+2 2Y+2

После приведения подобных членов линейная комбинация функций Чебышева (3) при-

Y

нимает вид Л^) = Bvuv. Эта функция является рабочей формой представления моде-

v =0

ли при ее непосредственном использовании для расчетов; число Y может рассматриваться как размерность модели.

С учетом соотношения (2) формула (1) для функции Грина упрощается:

1 т В

ОД «Т^У , . (4)

4пес \А2 + (ит)2 ' '

Соотношение (4) должно выполняться с необходимой для практического использования точностью.

3. Точность вычисления функции Грина

При решении задач, включающих решение уравнений Пуассона и Лапласа, традиционно ориентируются на погрешности выполнения законов сохранения и расчета полей менее 1% [7]. Оценим точность представления модели, основываясь на точности вычисления функции Грина по формуле (4). Для этого представим выражение (1) в виде

^ =уП- (/! + /2),

где

в те

/1 = I /о(АГ)Ф(А) ¿А; /2 = I /о(АГ) Ф(А) ¿А;

о в

в — произвольный предел интегрирования, выбираемый из условия

Ф(в) « Ф(го). (5)

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

/2 « у Jo(Ar) dA = Ф(го) ■ в(вг), (6)

в

где аналитическое выражение функции в(вг) имеет вид

в(вг) = 1 - вг Jo(er) + [Jo(er)Hi(er) - Ji(er)Ho(er)]. (7)

Здесь JK — функция Бесселя первого рода к-го порядка; HK — функция Струве к-го порядка; к = 1, 2 [8].

Точная верхняя граница погрешности вычисления /2 по формуле (6) sup А/ = [Ф(то) — Ф(в)]в(вО, что следует из теоремы о среднем [6] при монотонном характере изменения функции Ф(А).

Функция в(вО имеет колебательный характер (рис. 1). Первые нули 9К (к =1,16) этой функции представлены в табл. 1. При вычислении нулей с помощью выражения (7) функции Бесселя рассчитывались по интегральной формуле Бесселя [9], функции Струве — по интегральной формуле Пуассона [9]. При численном интегрировании использовалась

те

Таблица!. Нули ©-функции

к вк

1 1.108364661

2 4.062644472

3 7.151557848

4 10.269381067

5 13.397636191

6 16.530742102

7 19.666476075

8 22.803787326

9 25.942117194

10 29.081142112

11 32.220662040

12 35.360546791

13 38.500708053

14 41.641083834

15 44.781629355

16 47.922311469

Рис. 1. ©-функция.

квадратурная формула Гаусса для трех узлов [6], причем для снижения влияния методической погрешности интервал интегрирования был разбит на 104 шагов. Для снижения влияния погрешности округления все вычисления проводились с учетом 32 десятичных знаков мантиссы каждого операнда.

Соотношение (5) с заданной относительной погрешностью не более обеспечивается любым пределом интегрирования вк £ [в, Для выбора предела интегрирования интеграла 12 необходимо решить относительно в уравнение 1 —Ф(в)/Ф(то) = Определяемый по табл. 1 для наибольшего расстояния1 гтах ближайший больший нуль в-функции 9к даст произведение вттГтах, обладающее следующим свойством: для всех расстояний г < гтах может быть взят один и тот же нуль в-функции так как он соответствует значению предела вк = $к/г > втп, а вк заведомо обеспечивает соотношение (5). Для всех г < гтах при выбранном 9к < вттГтах предел вк всегда обеспечивает равенство 12 = 0 с требуемой точностью и при расчете Д может быть принята величина вк = $к/г. Тогда

вк

/ /о(Лг)Ф(Л) ¿Л. (8)

о

Вычисление интеграла с конечными пределами в (8) не представляет принципиальных трудностей при численном интегрировании. Значение функции Грина может быть рассчитано по формуле (8) практически с любой требуемой точностью, что позволяет оценить погрешность формулы (4), обусловленную аппроксимацией функции Ф(Л), т.е. неточностью представления модели.

Весьма существенное влияние на точность аппроксимации функции Ф(Л) оказывает выбор масштабирующего множителя т. В работе [4] предлагается принять множитель т

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

равным толщине одного из слоев (например, толщине подложки). В отсутствие формального критерия для выбора нужного слоя задачу определения т приходится решать методом полного перебора, что не всегда приемлемо из-за возникающих при вычислениях проблем: переполнение разрядной сетки, потеря точности и др. Более того, найденное таким способом значение т не позволяет обеспечить максимальную точность аппроксимации, которая достижима для заданного Т. Поэтому в качестве рабочего значения масштабирующего множителя предлагается принять величину то, которая минимизирует максимальную величину абсолютной погрешности аппроксимации функции Ф(А). Значение т = т0 назовем характеристическим размером конструкции.

4. Пример расчета

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

Ф(А)

где и = ехр(-Ато);

2 £1

— £1

1 1

1 + и2(^-Я)/то 1 + и2(П-Н)/то

(1+ и2л/то )[1 — и2(^-Я )/то ]_ (1 — и2Н/то )[1 + и2(^-я)/то]!

£1 и £2 — относительные диэлектрические проницаемости верхнего и нижнего слоев соответственно; Н — толщина нижнего слоя (толщина подложки); Н — расстояние между проводящими плоскостями (внутренняя высота корпуса).

В табл. 2 сведены результаты исследования точности математической модели конструкции в зависимости от ее размерности при т = Н и т = т0 для часто встречаемых на практике значений параметров модели: £1 = 1, £2 = 8.5, Н = 1 мм, Н = 3 мм. Минимум максимальной абсолютной погрешности аппроксимации математической модели Д0 и минимум максимальной относительной погрешности соответствуют величине т = т0,

Таблица2. Зависимость максимальной погрешности аппроксимации функции Ф(А) от размерности математической модели конструкции Т

т Д0 ¿0, % А 5, %

1 7.81-10-3 4.11 3.1410-2 18

2 1.2710-3 4.76 1.12-10-2 5.79

3 —1.6910-4 —1.21 9.5710-4 3.64

4 1.2910-5 1.53-10-1 —3.84-10-4 —1.8410-1

5 —3.93-10-6 —2.95-10-3 —5.910-5 —510-1

6 —1.5310-6 —1.0210-3 1.35-10-5 6.4-10-3

7 4.07-10-7 3.49-10-4 3.06-10-6 3.97-10-3

8 8.83-10-8 6.7-10-5 —6.510-7 —3.110-4

9 1.7310-8 1.2-10-5 —1.9410-7 1.48-10-4

10 3.9-10-9 1.9-10-6 2.9-10-8 1.38-10-5

11 —4.3-10-10 —3.28-10-7 —1.310-8 —6.9-10-6

Зв, % 2

О

-2

-4

-б

-8

-10

О 1 2

Рис. 2. Зависимость погрешности вычисления функции Грина от расстояния г при фиксированной размерности математической модели конструкции Т.

максимальные абсолютная А и относительная 8 погрешности рассчитывались при т = Н. Из таблицы видно, что при т = то точность аппроксимации функции Ф(А) значительно выше.

При Т = 2 получаем т0/Н = 1.842468, при Т = 3 — 1.970179. Для всех Т > 3 отношение т0/Н ~ 2 (с абсолютной погрешностью не более 6.5• 10-4). При увеличении е2 максимальная погрешность аппроксимации функции Ф(А) несколько возрастает, при увеличении Н — почти не меняется.

Соотношение (5) с относительной погрешностью 8ф = 10-2 % выполняется при в = 4896.11. При г = т0 получим вг = 9.64622. Ближайший больший корень в-функции в4 = 10.269381, при г = 2т0 — 96 = 16.530742. График на рис. 2 показывает погрешность вычисления функции Грина в плоскости, где расположен заряд, в зависимости от расстояния от заряда до точки, где вычисляется поле, при различных значениях Т. Все вычисления проводились с учетом 32 десятичных знаков мантиссы каждого операнда.

Выводы

1. Погрешность расчета функции Грина при традиционном методе вычислений с использованием формулы Вебера — Липшица существенно зависит от расстояния г.

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

3. Погрешность расчета функции Грина с использованием предлагаемого метода не зависит ни от расстояния, ни от параметров математической модели конструкции.

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

эвристическим обоснованием указанной размерности.

Список литературы

[1] Ollendorf F. Erdström. Berlin: Springer, 1928. 260 s.

[2] Плновский В., Филипс М. Классическая электродинамика. М.: Физматгиз, 1963. 432 с.

[3] Семенцов В.И. Расчет емкостей плоских проводников в слоистых средах // Радиотехника. 1973. Т. 28, № 10. С. 84-90.

[4] Тихонов А.Н., Самарский А.А. Уравнения математической физики. М.: Наука, 1977. 736 с.

[5] Конников И.А. Расчет емкостей прямоугольных пленочных проводников с произвольным коэффициентом формы // Судостроение. 1980. № 8. C. 32-33.

[6] Корн А.Г., Корн Т.М. Справочник по математике для научных работников и инженеров. СПб.: Лань, 2003. 831 с.

[7] Кукшева Э.А., Снытников В.Н. Параллельная реализация фундаментального решения уравнения Пуассона // Вычисл. технологии. 2005. Т. 10, № 4. С. 63-71.

[8] Градштейн И.С., Рыжик И.М. Таблицы интегралов, сумм, рядов и произведений. М.: Физматгиз, 1971. 1108 с.

[9] Handbook of mathematical functions / Ed. by M. Abramowitz, I.A. Stegun. N.Y., 1964. [Русский перевод: Справочник по специальным функциям / Под ред. М. Абрамовица, И. Стиган. М.: Наука, 1979. 400 с.]

Поступила в редакцию 4 апреля 2006 г., в переработанном виде —19 июня 2006 г.