ДИФФЕРЕНЦИАЛЬНЫЕ УРАВНЕНИЯ, 2022, том 58, № 8, с.1148-1157
ЧИСЛЕННЫЕ МЕТОДЫ
УДК 519.63
ВЫЧИСЛЕНИЕ ИНДУКТИВНОСТИ
НОРМАЛЬНЫХ ПРОВОДНИКОВ И СВЕРХПРОВОДНИКОВ
© 2022 г. М. М. Хапаев, М. Ю. Куприянов
Рассматривается задача вычисления матриц погонных индуктивности и сопротивления,
а также распределения тока и магнитного потенциала в протяжённой однородной линии
передачи, структура которой может состоять из сверхпроводников и нормальных провод-
ников произвольного сечения, характеризуемых некоторой комплексной проводимостью.
Даётся математическая постановка задачи и формулируется вычислительный алгоритм
граничных элементов, основанный на редукции задачи к системе граничных интеграль-
ных уравнений. Приводятся вычислительные альтернативы данному подходу и результа-
ты расчётов по разработанной программе TLZ (Transmission Line Z), а также предлагается
простой способ визуализации решений метода граничных элементов с помощью графиче-
ского пакета Gnuplot.
DOI: 10.31857/S0374064122080143, EDN: CHESMG
Введение. В работе рассматривается задача математического моделирования распределе-
ния электрического тока и электромагнитных полей в сечении двумерных проводящих струк-
тур, которые являются протяжёнными и однородными вдоль своей длины, так что проблема
сводится к решению двумерной краевой задачи в плоскости поперечного сечения. Образующие
структуру проводники могут быть нормальными металлами, сверхпроводниками и любыми
другими материалами, характеризуемыми некоторой комплексной проводимостью.
Интерес к таким задачам обусловлен развитием технологий сверхпроводниковой микро-
электроники [1-4]. Целью данной работы является разработка вычислительного метода и до-
статочно универсальной программы для вычисления матриц погонных частотно-зависимых
индуктивности и сопротивления произвольных структур, характеризуемых комплексной про-
водимостью.
Рассматриваемая задача может быть решена с помощью метода конечных элементов или
конечных разностей (см. [5, 6]) путём сведения к объёмному интегральному уравнению [7-9]
или c применением метода граничных элементов [7, 10, 11]. Ни один из этих подходов не явля-
ется стандартным, так как изучаемая задача, как это будет показано, является нелокальной.
Однако прямой метод граничных элементов обладает рядом преимуществ, например, точным
учётом скин-эффекта при расчёте пространственного распределения плотности тока в попе-
речном сечении изучаемых структур.
Ранее метод граничных элементов был применён для вычисления индуктивности сверх-
проводниковых линий передач без учёта диссипативных потерь [7, 11]. В данной статье метод
граничных элементов совершенствуется для вычисления комплексных частотно-зависимых
матричных индуктивных и резистивных составных частей комплексного импеданса исследу-
емых объектов. Исследовано совершенствование этого алгоритма с помощью использования
частичной заполненности матриц сеточных уравнений. Показано, как можно визуализировать
поля и токи любого метода граничных элементов с помощью технологий известной программы
Gnuplot [12].
1. Общая постановка физической задачи. Математическая модель задачи основана на
уравнениях Максвелла для комплексных амплитуд (j =
√-1) электрического поля
E, маг-
нитного поля
H , магнитной индукцииB,
B= μrμ0 H, векторного потенциала
A,
B= ∇×A
в калибровке ∇·A = 0, ток
J изаконаОм
J = σ(ω)E, где ω является угловой частотой, а
σ(ω) - комплексной проводимостью, μr - относительная магнитная проницаемость, которую
будем считать кусочно-постоянной и равной единице вне проводников. Для упрощения будем
пренебрегать током смещения.
1148
ВЫЧИСЛЕНИЕ ИНДУКТИВНОСТИ
1149
Пусть κ2 = -jωμrμ0σ(ω), тогда для векторного потенциала
A в проводнике получим
уравнение
2
∇×∇×A-κ2A=-κ
∇V, ΔV = 0,
(1)
jω
где V - потенциал электрического поля. В окружающей проводники среде уравнение для
векторного потенциала имеет вид
∇ × ∇ × A = 0.
Задача должна быть дополнена условиями на границах Γ сред с различной относительной
магнитной проницаемостью μr,
n - нормаль к границе:
[n ×A] = 0,
[n × (μ-1r∇ ×A)] = 0.
(2)
При описании сверхпроводников будем использовать двухжидкостную модель, в которой
проводимость задаётся выражением
1
σ(ω) = σn +
,
jωμ0λ2
где λ - лондоновская глубина проникновения [13], σn - нормальная проводимость [13]. Более
сложные модели проводимости в сверхпроводниках приведены, например, в [5], [14].
Для нормальных проводников известна модель Друде и её обобщения, согласно которой
проводимость определяется по формуле
σn
σ(ω) =
,
1-jωμ0τ
где τ - характерное время рассеяния электрона на примесях.
2. Постановка двумерной задачи для векторного потенциала. Получим из уравне-
ния (1) двумерную постановку задачи.
Будем считать, что имеется N +1 протяжённый вдоль оси z и однородный в сечении (x, y)
проводник. Проводник 0 служит для замыкания всех токовых петель и, таким образом, несет
возвратный ток. Каждый проводник не обязательно является односвязным. Каждая компонен-
та связности характеризуется своей постоянной комплексной проводимостью и соответственно
параметром k, а также относительной магнитной проницаемостью μm.
Пусть векторный потенциал
A и ток Jm имеют только z-компоненту:
J = {0,0,Jm(x,y)},
A= {0,0,A(x,y)}.
(3)
Пусть Sm - сечение каждого проводника с границей Γm, нормаль n к которой будем
считать внешней, r = (x, y) - радиус-вектор в плоскости поперечного сечения структуры,
Am(r) - z-компонента векторного потенциала при r ∈ Sm и A(r) в остальном пространстве,
Im - соответствующий полный ток, определяемый формулой
∫∫
Im = Jm(r)dxdy, m = 0,N.
Sm
Введём функцию ϕ = Vm/(jω), которая для сверхпроводников имеет смысл функции
фазы их параметров порядка [13]. Из (1) и (3) следует, что ∇xyϕm(x, y, z) = 0. Справедливо
следующее
Утверждение 1. Для каждого проводника m (m = 0, N ) имеет место равенство
ϕm(x, y, z) = ϕm(z).
В двумерном случае функция ϕm(x, y, z) не зависит от переменных (x, y). Это обстоя-
тельство позволяет значительно упростить постановку задач работ [5, 10, 15-19].
ДИФФЕРЕНЦИАЛЬНЫЕ УРАВНЕНИЯ том 58
№8
2022
1150
ХАПАЕВ, КУПРИЯНОВ
Далее положим
ϕm(z) = zΦm, Φm ≡ const,
(4)
что соответствует однородному распределению тока и электрического потенциала по длине
проводников и условию гармоничности потенциала ΔV = 0.
Можно определить зависимость от переменной z иначе (см. [5, 15]):
ϕm(z) = e-γzΦm, Φm ≡ const.
(5)
Комплексный параметр разделения переменных γ согласно работе [15] определим позднее.
В случае зависимости (4) постановка задачи относительно компоненты векторного потен-
циала A(r) имеет вид
ΔAm(r) + κ2mAm(r) = κ2mΦm, r ∈ Sm, m = 0,N,
(6)
⋃
ΔA(r) = 0, r ∈ R2 \
Sm,
(7)
m=0
1
∂Am(r)
∂A(r)
Am(r) = A(r),
=
,
r∈Γm,
(8)
μm
∂n
∂n
∫
1
∂Am(r)
ds = μ0Im, m = 0, N ,
(9)
μm
∂n
Γm
∑
Im = 0.
(10)
m=0
Условия сопряжения (8) следуют из (2). В равенстве (10) токи в проводниках Im, m = 1, N ,
замыкаются с помощью тока I0. Следствием баланса токов (10) является следующее
Утверждение 2. В задаче (6)-(10) справедливо предельное равенство
lim A(r) = 0.
(11)
r→∞
С помощью потенциала Am можно вычислить плотность тока в проводниках следующим
образом:
μ0Jm(r) = κ2m(Am(r) - Φm).
(12)
Исключив в (6) константу Φm с помощью формулы (12) (см. [20]), запишем это уравнение
в виде
∫∫
κ2m
κ2m
ΔAm(r) + κ2mAm(r) -
Am(r)ds +
Im = 0, r ∈ Sm.
(13)
|Sm|
|Sm|
Sm
В поставленной задаче (6)-(10) полные токи Im удобно считать заданными, а величины
Φm искомыми, в результате чего легко устанавливается связь между ними.
Заметим, что если выполняется зависимость (5), то постановка задачи остаётся такой же,
но для функций A(r)/γ и Am(r)/γ.
Величины Φm имеют смысл погонного потока магнитного поля или флюксоида в сверх-
проводимости. С помощью Φm нетрудно вычислить индуктивность и сопротивление, которые
в случае нескольких проводников имеют форму матриц.
3. Матрица импедансов. Традиционно индуктивность и сопротивление находятся с по-
мощью вычисления полей и токов [9, 15] из энергетических соображений. Между тем при
полной постановке задачи (6)-(10) с заданными токами Im и искомыми величинами Φm мат-
рица импедансов может быть определена более простым и точным способом.
ДИФФЕРЕНЦИАЛЬНЫЕ УРАВНЕНИЯ том 58
№8
2022
ВЫЧИСЛЕНИЕ ИНДУКТИВНОСТИ
1151
Так как задача (6)-(10) линейна, то существует такая комплексная матрица M размерно-
сти N × N, что
∑
Φi - Φ0 = MijIj.
j=1
Матрица импедансов Z(ω) = jωM(ω) имеет следующий вид:
Z(ω) = jωM(ω) = R(ω) + jωL(ω),
где R(ω) и L(ω) являются действительными N × N -матрицами погонных сопротивлений и
индуктивностей. Если размеры структуры в сечении приведены в микронах, то физической
размерностью этих матриц являются ом/микрон и пикогенри/микрон.
В работе [11] матрица индуктивностей для сверхпроводников получена с помощью функ-
ционала полной энергии.
Основной задачей вычислений настоящей работы является матрица комплексных импе-
дансов Z(ω), которая может быть использована как составная часть модели телеграфных
уравнений протяжённой линии передачи (см. [15]) или иным образом.
4. Постановка задачи для вычислительного алгоритма. Задача (6)-(10) обладает
рядом особенностей, которые должны быть учтены при численном её решении. Особенностью
постановки задачи являются неизвестные величины Φm, которые должны быть определены с
помощью дополнительных соотношений (9). Другая особенность - наличие малого параметра.
Этим параметром в случае сверхпроводников является лондоновская глубина проникнове-
ния λ, которая может принимать значение 0.1 микрон, в то время как толщина и ширина
проводников могут быть порядка единиц и десятков микрон, что приводит к существенной
концентрации тока в окрестности границ на расстояниях порядка λ.
Аналогичная проблема может возникнуть в случае нормальных металлов на больших ча-
стотах ω, где также возникает заметная концентрация тока вблизи границ проводников. Вы-
числительный алгоритм должен достаточно точно моделировать эти особенности.
Задача (6)-(10) может быть численно решена с помощью метода конечных элементов [5]
или конечных разностей [6]. Для этого уравнения (6) следует записать в более удобном для
решения виде - в виде интегро-дифференциальных уравнений (13).
Необходимость в таком подходе может возникнуть в случае пространственной неоднород-
ности проводимости σ(ω, r). Такой подход требует построения алгоритма разрешения зависи-
мости Φm от Im, что приводит либо к итерационному процессу [5, 6], либо требует окаймления
матрицы сеточных уравнений строками и столбцами для аппроксимации соотношений с Φm
и Im в случае алгебраического подхода к сборке системной матрицы. Что касается сеток,
то проблема концентрирования тока приводит к их сильному сгущению, которое показано в
статье [21]. Существует также проблема выбора достаточно большого вмещающего окна для
решения задачи во всем пространстве.
Другим простым и достаточно очевидным вычислительным методом является сведение
задачи к уравнению по объёму проводника [9, 15], которое в случае задачи (6)-(10) приведёт
к объёмному интегральному уравнению по сечениям проводников
∫∫
(
)
∑
1
1
1
Jm(r0) +
ln
Jl(r)dxdy - Φm = 0,
(14)
κ2m
2π
|r - r0|
l=0
Sl
∫∫
Jm(r)dxdy = Im, m = 0,N.
(15)
Sm
В отличие от работы [15] данная постановка явно включает константы Φm и соотношения
полного тока (15).
Численное решение уравнения (14) требует дискретизации сечений проводников и при-
менения какого-либо проекционного метода: коллокации или Галёркина. Проблемой в этом
ДИФФЕРЕНЦИАЛЬНЫЕ УРАВНЕНИЯ том 58
№8
2022
1152
ХАПАЕВ, КУПРИЯНОВ
случае является очень большая размерность матрицы системы сеточных уравнений (см. [7]),
вызванная двумерной дискретизацией областей Sm. Вместе с тем проблема концентрирова-
ния тока требует применения достаточно мелких сеток для точной аппроксимации решения
вблизи поверхности проводников. Имеется также трудность, связанная с недостаточной чис-
ленной устойчивостью такого алгоритма, так как коэффициент при внеинтегральном члене
обычно мал и задача близка к классической некорректной задаче.
Другим известным подходом является алгоритм, в котором объёмное интегральное урав-
нение (14) преобразуется к интегральному уравнению по границе проводников со сложным
нестандартным ядром [16] и новой плотностью, относительно которой записывается уравне-
ние по поверхности проводников. Согласно этому подходу внутри каждого проводника для
плотности тока используется представление типа потенциала простого слоя, а z-компонента
электрического поля удовлетворяет уравнению Гельмгольца ΔEz,m + κ2mEz,m = 0, решение
которого представимо в виде потенциала простого слоя с плотностью s(q) и ядром Mm(p, q):
∮
Jm(p)
jπ
Ez,m(p) =
= Mm(p,q)s(q)dlq, Mm(r,r0) = -
H(2)0(κm|r - r0|).
σ(ω)
2
Γm
Подставив выражение Jm(r) = σ(ω)Ez,m(r) в (14), получим интегральное уравнение с
ядром
∫∫
1
Gk(r0, p) =
ln
Mk(r,p)dsr.
(16)
|r - r0|
Sk
Ядро (16) представляет собой интеграл по сечению проводников, который в работе [16] пред-
лагается вычислять с помощью квадратур и сетки, введённой в Sm, что значительно увели-
чивает время на формирование матрицы системы сеточных уравнений, полученных тем или
иным способом.
Наконец, для задачи (6)-(10) может быть сформулирован прямой метод граничных эле-
ментов. Для сверхпроводников этот подход ранее был применён в [7, 11]. В данной работе он
перенесён на случай проводников с комплексной проводимостью. Известны и другие реализа-
ции прямого метода граничных элементов для случая нормальных проводников [10, 18, 20],
в отличие от которых в настоящей работе явно используются величины Φm и Im и простой
метод вычисления матрицы импедансов с помощью потоков Φm.
5. Интегральные уравнения прямого метода граничных элементов. Прямой метод
граничных элементов состоит в применении третьей формулы Грина для редукции дифферен-
циальной задачи (6)-(10) к интегральным уравнениям по границе проводников Γm.
Интегральные уравнения прямого метода граничных элементов относительно искомых
функций Am(r) и ∂Am(r)/∂n, где n(r) - внешняя нормаль к Sm, а также искомых величин
Φm, имеют вид
∫
(
)
∑
∂G(r, r0)
1 ∂Al(r)
πAm(r0) +
Al(r)
-
G(r, r0) dlr = 0,
(17)
∂n
μl
∂n
l=0
Γl
∫
(
)
∂Mm(r,r0)
∂Am(r)
πAm(r0) -
Am(r)
-
Mm(r,r0) dlr +
∂n
∂n
Γm
(
∫
)
∂Mm(r,r0)
+ π-
dlr Φm = 0, r0 ∈ Γm,
(18)
∂n
Γm
∫
1
∂Am(r)
dlr = μ0Im, m = 0, N,
(19)
μm
∂n
Γm
ДИФФЕРЕНЦИАЛЬНЫЕ УРАВНЕНИЯ том 58
№8
2022
ВЫЧИСЛЕНИЕ ИНДУКТИВНОСТИ
1153
1
jπ
G(r, r0) = ln
,
Mm(r,r0) = -
H(2)0(κm|r - r0|),
(20)
|r - r0|
2
где r0 - точка на Γm, отличная от угловой.
В этих уравнениях G(r, r0) - фундаментальное решение уравнения Лапласа (7), Mm(r, r0) -
фундаментальное решение уравнения Гельмгольца (6) (Im (κm) < 0).
Для сверхпроводника при ω = 0 имеем Mm(r, r0) = K0(|r - r0|/λ), где K0(x) - модифи-
цированная функция Бесселя нулевого порядка.
Уравнение (17) является третьей формулой Грина для неограниченной области (см. [22]), в
которой учтено условие (11), следующее из баланса полных токов (10), а также равенства (8).
Уравнение (18) является третьей формулой Грина для сечения проводника Sm, записанной
для функции u(r) = Am(r) - Φm следующим образом:
∫
(
)
∂Mm(r,r0)
∂Am(r)
α(r0)(Am(r0) - Φm) -
(Am(r) - Φm)
-
Mm(r,r0) dsr = 0,
∂n
∂n
Γm
где α(r) (r ∈ Γm) - угол видимости кривой Γm из точки r.
Окончательно приходим к уравнениям (17)-(20), в которых заданы полные токи Im, а
величины Φm, наряду с функциями Am(r) и ∂Am(r)/∂n, на границах проводников Γm яв-
ляются искомыми.
Задав Im = δkm (m, k = 1, N ) и решив поставленную задачу (17)-(20), можно построчно
вычислить матрицу импедансов Z = R + jωL.
6. Вычислительный метод для системы граничных интегральных уравнений.
В вычислительном алгоритме данной работы и программе граница каждого из проводников
имеет следующую составную структуру:
⋃
Γk = γi,
i=1
где каждая из элементарных кривых γi может быть отрезком прямой, дугой окружности,
окружностью либо кривой второго порядка, проходящей через три заданные точки. Все эти
элементарные кривые не имеют общих точек, кроме концов.
Далее, на каждой из кривых γi вводится сетка, которая может быть как равномерной,
так и неравномерной, причём угловые точки, являющиеся концами γi, узлами сетки не явля-
ются. На этой сетке аппроксимируются уравнения (17)-(20) с помощью метода коллокации.
Интегралы в соотношениях (19) для полного тока аппроксимируются с помощью квадрату-
ры трапеций. Искомыми величинами являются функции Am(r) и (∂Am(r)/∂n)/μm) в узлах
сетки, а также константы Φm.
Более подробно аппроксимация интегральных уравнений описана в статьях [7, 11]. В на-
стоящей работе используется тот же подход.
Вычислительный метод приводит к системе линейных сеточных уравнений Ku = f, в пра-
вой части которой оказываются заданы полные токи Im. Матрица K не имеет специальной
структуры, однако обладает некоторой разреженностью, которая может быть использована
при численном решении. Структура этой разреженности хорошо известна и обусловлена тем,
что векторный потенциал A(r) в каждом проводнике не связан непосредственно с потенциа-
лом в других проводниках. Все проводники в единую систему уравнений объединяет потенциал
во вмещающей среде.
Основной целью вычислений является матрица импедансов, но также важной задачей яв-
ляется визуализация векторного потенциала и распределения тока. Для этого можно исполь-
зовать известную программу Gnuplot [12]. Обычно эта программа используется для построе-
ния графиков, в том числе достаточно сложных [23]. В [24] показано как можно применить
Gnuplot для визуализации конечно-элементных решений, использовав вычислительную сет-
ку и правильную раскраску ячеек в зависимости от значений решения. В рассматриваемом
в работе случае сетка в области решения не строится, что является преимуществом метода
10 ДИФФЕРЕНЦИАЛЬНЫЕ УРАВНЕНИЯ том 58
№8
2022
1154
ХАПАЕВ, КУПРИЯНОВ
граничных интегральных уравнений. Тем не менее можно достаточно просто построить визу-
ализацию решения с помощью алгоритма dgrid3d интерполяции нерегулярно распределённых
данных, реализованного в Gnuplot.
В методах граничных интегральных уравнений легко реализуется вычисление решения в
заданной точке пространства. Векторный потенциал A(x, y) в случае вакуума вычисляется
по формуле
∫
)
∑
1
( 1 ∂Al(r)
∂G(r, r0)
A(r0) =
G(r, r0) - Al(r)
dlr.
(21)
2π
μl
∂n
∂n
l=0
Γl
Для каждого проводника справедливо представление
∫
(
)
1
∂Mm(r,r0)
∂Am(r)
Am(r0) =
(Am(r) - Φm)
-
Mm(r,r0) dlr + Φm, r0 ∈ Sm.
(22)
2π
∂n
∂n
Γm
Интегралы в формулах (21), (22) аппроксимируются численно, аналогично аппроксимации
уравнений (17)-(20). С помощью этих выражений можно вычислить A(r), Am(r) и Jm(r) по
формуле (12) на различных сетках, пропустив точки вблизи или на границах Γm. В этом
случае ядра интегральных представлений решения (21), (22) не сингулярны и численно ин-
тегрируются с достаточной для визуализации точностью. Кроме того, при решении задачи
вычисляются Am(r) в определённых точках Γm. Все эти величины со своими координатами
образуют нерегулярно распределённые данные, по которым можно построить то или иное гра-
фическое представление результатов. Пример таких сеток в простом случае одного проводника
в пространстве приведён на рис. 1.
Рис. 1. Сетка с узлами трёх видов, используемых для визуализа-
ции. Узлы, отмеченные квадратами, используются для значений во
вмещающей среде и вычисляются с помощью формулы (21). Значе-
ния, отмеченные кружками, вычисляются при решении граничных
интегральных уравнений. Узлы, отмеченные треугольниками, со-
держат значения, вычисленные с помощью (22).
Таким образом, все данные, вычисленные на различных сетках, могут быть переданы и
нарисованы программой Gnuplot с помощью алгоритма интерполяции нерегулярных данных.
Этот подход для визуализации полей применим для любого метода граничных интеграль-
ных уравнений.
Описанный вычислительный подход был реализован в C++ программе TLZ (Transmission
Line Z), которая может быть скомпилирована для операционных систем Windows и Linux. Для
работы с разреженными матрицами использовался пакет Eigen [25].
7. Результаты расчётов. Для проверки точности алгоритма и программы были про-
ведены расчёты для нормальных проводников, имеющих коаксиальную структуру. В этом
случае известно аналитическое решение, которое воспроизводилось программой с хорошей
ДИФФЕРЕНЦИАЛЬНЫЕ УРАВНЕНИЯ том 58
№8
2022
ВЫЧИСЛЕНИЕ ИНДУКТИВНОСТИ
1155
точностью. Прочие результаты для нормальных проводников были опубликованы в [19] и по-
лучены несколько другим методом граничных интегральных уравнений, в котором условия
Φm = const учитывались не вполне точно. Рассматривались два симметрично расположенных
проводника с размерами сечений [0, 2]×[-0.6, -0.4] мм и [0, 2]×[0.4, 0.6] мм и с проводимостью
56 MСм/м. Результаты расчётов приведены в таблице.
Таблица. Значение параметров для двух широких тонких провод-
ников с сечениями [0, 2] × [-0.6, -0.4] мм и [0, 2] × [0.4, 0.6] мм и с
проводимостью 56 MСм/м
TLZ
Метод [20]
Частота, Гц
R, МОм/м L, нГн/м R, МОм/м L, нГн/м
103
89.1
357
87.6
357.1
106
181.7
335
175.3
333.9
109
5688
311
5790
310
Хотя были использованы различные вычислительные методы и подбора сетки не произво-
дилось, результаты совпали с относительной точностью 3.6 % для сопротивлений R и менее
1% для индуктивностей L.
8. Простая сверхпроводниковая структура. Рассмотрим известный пример (см. ра-
боты [4, 26]) структуры из двух сверхпроводящих прямоугольных проводников, нижний из
которых несет возвратный ток (рис. 2). В [26] в одной из задач проводники имеют размеры
[-4, 4] × [0.2, 0.5] мкм (нижний) и [-1.25, 1.25] × [0.6775, 0.8975] мкм (верхний). Для верхнего
проводника лондоновская глубина проникновения λ = 0.137 мкм, для нижнего - 0.086 мкм.
Результат измерения индуктивности в работе [26]: 0.1599 пГн/мкм, результат данной статьи -
0.1593 пГн при шаге дискретизации интегральных уравнений 0.25 мкм.
Рис. 2. Распределение действительной части z-компоненты векторного по-
тенциала A(r) (сверху) и тока (снизу) - чем темнее область, тем больше
значение потенциала. Размеры представлены в мкм. Для визуализации ис-
пользовалась программа Gnuplot c опцией dgrid3d (интерполяция нерегу-
лярно распределённых данных).
В настоящей работе на рис. 2 представлена немного другая структура. Для лучшего пред-
ставления скинирования тока толщина верхнего проводника была увеличена, так что по вы-
соте проводник занимал область [0.6775, 1.25], λ = 0.15 мкм. Результаты вычислений были
экспортированы в файл и на рисунке представлена визуализация этих данных.
Для задачи рис. 2 шаг сетки дискретизации интегральных уравнений 0.025 мкм привёл
к матрице размерности 1862 × 1862. Заполненность этой плотной матрицы сеточных урав-
нений составила 57%. Стандартный алгоритм LU факторизации (метод Гаусса) с выбором
ведущего элемента потребовал для нее 15.2 с. Затем та же матрица сеточных уравнений была
представлена в компактном разреженном формате, переупорядочена и факторизована с по-
мощью пакета Eigen [25]. Операция разреженной факторизации длилась 1.7 с. Таким образом,
ДИФФЕРЕНЦИАЛЬНЫЕ УРАВНЕНИЯ том 58
№8
2022
10∗
1156
ХАПАЕВ, КУПРИЯНОВ
для сеточных уравнений прямого метода граничных элементов с помощью аппарата разре-
женных матриц можно добиться значительного ускорения решения.
Заключение. В данной работе представлена постановка задачи для векторного потенци-
ала в протяжённых линиях передачи, состоящих из проводников с комплексной проводимо-
стью, в частности, сверхпроводников. Особенностью предложенного метода решения задачи
является использование интегро-дифференциальных уравнений для векторного потенциала,
содержащих величины типа интегральных погонных потоков, с помощью которых легко вы-
числяется матрица погонных импедансов.
Показано как свести поставленную задачу к граничным интегральным уравнениям по гра-
нице проводящих структур и ускорить численное решение этих интегральных уравнений, ис-
пользовав их определённую разреженность.
Простая визуализация решения методом граничных элементов получена с помощью про-
граммы Gnuplot с опцией dgrid3d.
В дальнейшем планируется усовершенствовать алгоритм для случая составных проводни-
ков и неоднородной вмещающей среды.
Работа выполнена при финансовой поддержке Российского научного фонда (проект 20-12-
00130).
СПИСОК ЛИТЕРАТУРЫ
1. Krylov G., Kawa J., Friedman E.G. Design automation of superconductive digital circuits: a review
// IEEE Nanotechnology Magazine. 2021. V. 15. № 6. P. 54-67.
2. Tolpygo S.K. Advanced fabrication processes for superconductor electronics: current status and new
developments // IEEE Transactions on Applied Superconductivity. 2019. V. 29. № 5. P. 1-13.
3. Oates D.E., Tolpygo S.K., Bolkhovsky V. Submicron Nb microwave transmission lines and components
for single-flux-quantum and analog large-scale superconducting integrated circuits // IEEE Transactions
on Applied Superconductivity. 2017. V 27. № 4. P. 1-5.
4. Tolpygo S.K., Golden E.B., Weir T., Bolkhovsky V. Inductance of superconductor integrated circuit
features with sizes down to 120 nm // Superconductor Science and Technology. V. 34. P. 085005.
5. Roux P.L., Jackman K., Delport J.A., Fourie C.J. Modeling of superconducting passive transmission
lines // IEEE Transactions on Applied Superconductivity. 2019. V. 29. № 5. P. 1-5.
6. Alsop L.E., Goodman A.S., Gustavson F.G., Miranker W.L. A numerical solution of a model for a
superconductor field problem // J. of Comput. Phys. 1979. V. 31. № 2. P. 216-239.
7. Хапаев М.М. Метод граничных интегральных уравнений для одной модели электрического тока в
сверхпроводниках // Журн. вычислит. математики и мат. физики. 1998. Т. 38. № 1. С. 115-121.
8. Sheen D.M., Ali S.M., Oates D.E., Withers R.S., Kong J.A. Current distribution, resistance, and
inductance for superconducting strip transmission lines // IEEE Transactions on Appiled Superconduc-
tivity. 1991. V. 1. № 2. P. 108-115.
9. Tsuk M.J., Kong J.A. A hybrid method for the calculation of the resistance and inductance of
transmission lines with arbitrary cross sections // IEEE Transactions on Microwave Theory and
Techniques. 1991. V. 39. № 8. P. 1338-1347.
10. Ruey-Beei Wu, Jin-Chum Yang. Boundary integral equation formulation of skin effect problems in
multiconductor transmission lines // IEEE Trans. on Magn. 1989. V. 25. № 4. P. 3013-3015.
11. Khapaev M.M. Extraction of inductances of a multi-superconductor transmission line // Supercond. Sci.
Technol. 1996. V. 9. P. 729.
13. Orlando T.P., Delin K.A. Foundations of Applied Superconductivity. Austin, 1991.
14. Zimmermann W., Brandt E., Bauer M., Seider E., Genzel L. Optical conductivity of BCS superconduc-
tors with arbitrary purity // Physica C: Superconductivity. 1991. V. 183. № 1-3. P. 99-104.
15. Plaza G., Marques R., Medina F. Quasi-TM MoL/MoM approach for computing the transmission-line
parameters of lossy lines // IEEE Transactions on Microwave Theory and Techniques. 2006. V. 54. № 1.
P. 198-209.
16. Menshov A., Okhmatovski V. Novel surface integral equation formulation for accurate broadband RL
extraction in transmission lines of arbitrary cross-section // IEEE/MTT-S Int. Microwave Symposium
Digest. 2012. P. 1-3.
ДИФФЕРЕНЦИАЛЬНЫЕ УРАВНЕНИЯ том 58
№8
2022
ВЫЧИСЛЕНИЕ ИНДУКТИВНОСТИ
1157
17. Coperich K.M., Morsey J., Okhmatovski V.I., Cangellaris A.C., Ruehli A.E. Systematic development
of transmission-line models for interconnects with frequency-dependent losses // IEEE Transactions on
Microwave Theory and Techniques. 2001. V. 49. № 10. P. 1677-1685.
18. Patel U.R., Triverio P. Skin effect modeling in conductors of arbitrary shape through a surface admittance
operator and the contour integral method // IEEE Transactions on Microwave Theory and Techniques.
2016. V. 64. № 9. P. 2708-2717.
19. Djordjevic A.R., Sarkar T.K., Rao S.M. Analysis of Finite Conductivity Cylindrical Conductors Excited
by Axially-Independent TM Electromagnetic Field // IEEE Trans. Microwave Theory and Techn. 1985.
V. 33. P. 960-966.
20. Konrad A. Integrodifferential finite element formulation of twodimensional steady-state skin effect
problems // IEEE Trans. on Magn. 1985. V. 18. P. 284-292.
21. Herbst H.F., Roux P.L., Jackman K., Fourie C.J. Improved transmission line parameter calculation
through tcad process modeling for superconductor integrated circuit interconnects // IEEE Transactions
on Applied Superconductivity. 2020. V. 30. № 7. P. 1-4.
22. Кошляков Н.С., Глинер Э.Б., Смирнов М.М. Уравнения в частных производных математической
физики. М., 1970.
24. Glanzel J., Unger R. High quality FEM-postprocessing and visualization using a gnuplot based toolchain
// Chemnitz Scientific Computing Preprints. 2014.
26. Chang W.H. The inductance of a superconducting strip transmission line // J. of Appl. Phys. 1979. V. 50.
№ 12. P. 8129-8134.
Московский государственный университет
Поступила в редакцию 05.03.2022 г.
имени М.В. Ломоносова,
После доработки 26.04.2022 г.
Научно-исследовательский институт ядерной физики
Принята к публикации 05.07.2022 г.
имени Д.В. Скобельцына (НИИЯФ МГУ), г. Москва
ДИФФЕРЕНЦИАЛЬНЫЕ УРАВНЕНИЯ том 58
№8
2022
