CC BY
28
Вычислительные технологии
Том 3, № 3, 1998
ОБ ОДНОМ ЧИСЛЕННОМ МЕТОДЕ НАХОЖДЕНИЯ СТАЦИОНАРНЫХ РЕШЕНИЙ ГИДРОДИНАМИЧЕСКОЙ МОДЕЛИ ПЕРЕНОСА НОСИТЕЛЕЙ ЗАРЯДА В ПОЛУПРОВОДНИКАХ*
А. М. Блохин, А. С. БУШМАНОВА Институт математики СО РАН, Новосибирск, Россия e-mail: blokhin@ma.th.nsc.ru, bushma.nova@math.nsc.ru
Numerical method for finding the stationary solutions of the hydrodynamic model equations for charge transport in semiconductors is described in the paper. With the help of this method based on generalized Green matrix the existence the approximate solution of the well-known ballistic diode problem is obtained.
1. Введение
В настоящее время для нахождения решений стационарных одномерных гидродинамических уравнений, описывающих физическую сущность явлений переноса носителей заряда в полупроводниках, используются самые различные вычислительные алгоритмы. Так, в работах [1-3] были применены специально сконструированные конечно-разностные методы для расчета стационарных режимов гидродинамических уравнений. В [4-6] для определения стационарных режимов одномерных гидродинамических уравнений используется принцип установления, при применении которого решение стационарной задачи находится как предел решения нестационарной задачи при £ ^ то .
В настоящей работе для расчета стационарного решения известной задачи о баллистическом диоде предлагается численный метод, основанный на использовании так называемой обобщенной матрицы Грина.
2. Предварительные сведения и постановка задачи
Рассмотрим газодинамическую модель переноса зарядов в полупроводниках, описываемую четырьмя уравнениями (более подробно эта модель обсуждается в [5], там же приводятся
* Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований, грант №96-01-01560.
© А. М. Блохин, А. С. Бушманова, 1998.
3
результаты численных исследований этой модели):
Ят + ^1} — 0,
41} + ^2) — R^s +
Дм
3
ДЕ - - Д
(ДЕ)т + . ) — Дм^ +-----------г—2— — 0,
т,,
¥ss — в (Д — Р),
(1)
(2)
(3)
(4)
где соответствующие потоки
3(1), 3(2), 3(3)
задаются формулами
3(1) = Дм,
3(2) = Дм2 + р,
3(3) = ДмП,
15 13
причем р = Д$ — давление; П = ^м2 + ^$5 Е = ^м2 + ^$ — полная энергия. Здесь неизвестные функции Д, $, м — соответственно плотность, температура и скорость электронного газа, ^ — электрический потенциал. Времена релаксации тр, тш считаются функциями полной энергии (см. [5]), р = р(в) — плотность легирования или естественная электронная концентрация (известная функция), в — константа.
Задача о баллистическом диоде рассматривается на интервале 0 < в < 1. Граничные условия имеют вид:
Д(0) = Д(1) = 1, ^(0) = А, ^(1) = А + В,
$(0) = 1, (5)
где А, В — константы, причем В > 0 — напряжение смещения (см. [4, 5]). В стационарном случае исходная система будет следующей:
(Дм).
0,
Дм
(Дм + — 0,
Дм
тр
Д( 1 м2 + 3 ^ — 3 Д 2 2)2
'Т'ІП
0
р
0
s
^ss — в (Д — Р)- (6)
Граничные условия здесь прежние. Приведем систему (6) к нормальному виду, т. е. разрешим ее относительно производных. При этом условие разрешимости записывается как
Дм(5# — 3м2) — 0.
Тогда (6) перепишется в следующем виде:
Q1 Яв = в(Д — р) Дв = ДФ, мв = —мФ,
1 м2 + 3$ — 3 2 2
$в = зЯ — 7------~------ ^Фм2,
5 5 тт м 5
(7)
где
Ф
1
5# — 3м2
_ 5м м2 + 3# — 3
зд — — +
т,„ м
В векторной форме задача о баллистическом диоде формулируется следующим обра-
зом:
2' — А2 + Т (в, 2),
MZ — MoZ (0) + Мі 2 (1) — В,
(8)
(9)
где
( <Р \
д Д
м #
А
( 0 1 0 0 0 \
0 0 в 0 0
0 0 0 0 0
0 0 0 0 0
0-000
/
( тЛ
Т2
, т — Тз —
Т4
V т5 /
0
—вр ДФ
—мФ 1 м2 + 3# — 3 5 Т,„ м
22 + -Фм2
5
1 0 0 0 0 0 0 0 0 0 \ / А \
0 0 1 0 0 0 0 0 0 0 1
Мо — 0 0 0 0 1 , Мі — 0 0 0 0 0 , В — 1
0 0 0 0 0 1 0 0 0 0 А + В
0 0 0 0 0 0 0 1 0 0 / V 1 /
Рассмотрим соответствующую “однородную” задачу
2' — А2, М2 — 0.
(10)
(11)
Очевидно, что фундаментальную матрицу для системы 2' = А% можно взять в виде
(
X (в)
в
1
0
00 2 2
V 5 5
2
в
1
в
0
в
01
00
00
10
00
\
р
2
в
2
в
5
Рассмотрим D = MX — матрицу, полученную подстановкой в краевое условие фундаментальной матрицы. Известно, что матрица Грина для однородной задачи существует и единственна тогда и только тогда, когда rank D = n, где n = dim X. В нашем случае
detD
1 0
0 0 a 0 0
0 0 1 1 0 0 00
1
11-00
00
2
1
в
00
т. е. rank D < n и обычную матрицу Грина построить нельзя, однако можно построить обобщенную матрицу Грина (см. [7]).
3. Обобщенная матрица Грина
Наряду с нелинейной задачей (8), (9) рассмотрим линейную неоднородную задачу
Zs = AZ + f (s), (12)
MZ = B. (13)
Известно,что общее решение системы линейных дифференциальных уравнений (12) имеет вид
Z (s) = X (s)c + Z(s), (14)
где Z(s) — частное решение неоднородной системы (12), которое мы определим как
Z(s) = / K(s,t)/(t) dt,
здесь К(в,і) — п х п-мерная матрица
К(в,і) — - X (в)Х-1 (і)зі^(в — і).
Чтобы вектор 2(в) был решением краевой задачи (12), (13), необходимо удовлетворить также краевым условиям. Подставляя (14) в (13), получим алгебраическую систему относительно вектора с Є Дп:
Dc + MZ = B,
(15)
где MZ = MoZ(0) + MiZ(1) = M / K(-, t)f (t) dt.
o
Пусть rank D = n1 < n; тогда, согласно [7], система (15) разрешима тогда и только тогда, когда
Pr*(B-MZ) = 0, r = n - Пі.
0
і
0
і
При этом общее решение системы (15) имеет вид
с = Есг + Р+(В-М),
где Рг* — матрица г х п, состоящая из линейно-независимых строк матрицы Р*, Е — фундаментальная п х г-мерная матрица решений однородной системы, соответствующей (15), сг — произвольный г-мерный вектор-столбец из Кг.
Поясним вышеприведенные обозначения. Пусть Т>* — п х п-мерная матрица, сопряженная к Р. Обозначим через N(Р) и N(Р*) нуль-пространства матриц (ядра операторов) V и Т>* = Т>т соответственно: N(Р) = кег V = {с : с € Яга, Рс = 0}, N(Р*) = кег Т>* = {^ : 5 € Дга, вР = 0}. Тогда Р и Р* — ортопроекторы, проектирующие Дга на нуль-пространства N(Р) и N(Р*). Р+ — единственная п х п-мерная матрица, псевдо-обратная (обобщенно-обратная) по Пенроузу к п х п-мерной матрице Р. Для нахождения ортопроекторов применяется следующая процедура. Пусть вектор-столбцы {/*},* = 1,..., г определяют базис нуль-пространства N(Р), а вектор-строки {$&}, к = 1,..., г — базис нуль-пространства N(Р*). По линейно независимым базисным векторам составляем матрицы
Грама а = {^ = (//= 1,...,г},7 = {7у = ^г.^' = 1,...,г}. Пусть а“1,7г^1 —
элементы матриц а 1,7 1. Тогда для ортопроекторов Р и Р* справедливы формулы
Г Г
Р = £ а/. Р* = £ ‘»зТ.
*,.7=1 г,^=1
Кроме того, здесь имеют место выражения, связывающие псевдообратную матрицу и ортопроекторы:
+
(Рт Р + Р) _1Рт
+
/„ - Р *, Р+Р = 1п - Р.
где /га — единичная матрица порядка п. В нашем случае
(
Р+
00
1 в -1 - 4
5
2
0
0
0
V
0
0
1
в
2
0
Р
00000 00000 00000 0 0 0 1 0 00000
00 00
0 - 20
0
1 - 4
в 2 0
Р*
0
0
0
0 0 0 0
1 0 0 -1
22
000 000 1
V
0 0 1
2 )
Подставляя найденный вектор с в (14), приходим к следующему выражению для общего решения задачи (12), (13):
2(в,сг) = Хг(в)сг + / К(з,£)/(¿) ^ - X(з)Р+М / К(^,£)/(¿) ^ + X(з)Р+В, (16)
0
1
1
о
о
где Хг (в) — п х г-мерная матрица, столбцы которой — полная система г линейно независимых решений однородной краевой задачи (10), (11): Хг(в) = X(в)Р.
Общее решение (16) краевой задачи (12), (13) состоит из суммы двух решений: общего решения Хг(в)От краевой задачи (10), (11) и частного решения неоднородной краевой задачи (12), (13). Выберем частное решение неоднородной краевой задачи (12), (13) из условия ортогональности на [0, 1] к любому решению краевой задачи (10), (11), т. е. из условия
1
J Хт(в)т2(в, ог) ¿в = 0. (17)
о
1
Обозначим через в = / Хг (в)тХг (в) ¿в — г х г-мерную невырожденную матрицу, через
о
1
втп = § Хг (в)тХ(в) ¿в — п х г-мерную. Из условия (17) единственным образом находим
о
вектор ог :
г 1 1 1
/ / Хг(в)тК(в,Ь) ¿в/(Ь) — 8гпТ>+М§ К(^,Ь)/(Ь)
о о о
-в-1 БтпЪ+В,
cr — —S
-і
при этом из множества решений Z(s,cr) вида (16) краевой задачи (12), (13) определяем единственное решение, удовлетворяющее условию ортогональности (17):
1
Zo(s) = J G(s, t)f (t) dt + [X(s) — Xr(s)S-1Sr„] D+B. (18)
0
Здесь приняты следующие обозначения. Если функционал M такой, что имеет место равенство
1 1 M j K(-,t)f (t) dt = JMK(-,t)f (t) dt,
0 0
то G(s,t) — n x n-мерная обобщенная матрица Грина краевой задачи (12), (13):
1
G(s, t) = G0(s, t) — Xr(s)S-1 J XГ(s)G0(s, t) dt,
0
G0(s,t) = K(s,t) — X (s)D+MK( ■ ,t).
n ( \ G01(s,t), s<t, . . \ G1(s,t), s<t, . ,
G°(M) = I Gm(s,t), s > t. С(аЛ) = ( G2(s,t), s > t. (19)
Итак, согласно [7], справедливо
Утверждение. Если rank D = n1, то краевая задача (10), (11) имеет r = n — n1 и только r линейно независимых решений. Неоднородная краевая задача (12), (13) разрешима тогда и только тогда, когда f (s) £ C[0,1] и B £ Rn удовлетворяют условию
P *
P r
B — MJ K(■ ,t)f (t) dt 0
0
і
и при этом имеет г-параметрическое семейство решений вида
2(в, сг) = Хг(в)сг + 2э(в),
где 20(в) определяется по (17).
Обобщенная матрица Грина С(в,Ь), определенная формулами (19), обладает следующими основными свойствами (проверка которых достаточно очевидна).
1. Каждый столбец матрицы при в = Ь является непрерывно дифференцируемым решением дифференциальной системы (10)
дС(в, Ь) дЬ
А(в)С(в,Ь), в = Ь.
2. В точке в = Ь матрица имеет разрыв по Ь первого рода
С(Ь + 0,Ь) — С(Ь — 0,Ь) = 1п.
3. Матрица удовлетворяет краевому условию
МС(■ ,Ь) = Р*МК(■ ,Ь), Р* : Кп ^ N(V*).
4. Каждый столбец матрицы ортогонален к любому решению краевой задачи (10), (11)
1
[ XТ(в)С(в,Ь) ¿в = 0.
Для рассматриваемой задачи
( 1 в — Ь 2
К(в,Ь) =
0
0
0
0
2
1
2
0
0
в — Ь 5
в
4
в
2
(в — Ь)2 (в — Ь)
1 2 0
& —Ь)!
00
00
00 1
2
0
0 1
2 )
sign(в — Ь),
о
МК( ■ ,Ь)
Ь
2
Ь
5
2
00
вЬ2
1
~2
То
1
2
0
0
0
0
0
0
1
2
00
0
0
0
0о1(в,Ь) =
—в
—1
0
в(Ь — 1) — (в + 2Ь2 — 4Ь + 1)
Ь — 1 — в(2в + 2Ь2 — 4Ь +1)
0
22
0
0
V----в
V 5
вв
5 в(Ь — 1) — 10(в + 2ь2 — 4Ь + 1)
00
00
00
1
2 0
00
002(в, Ь)
1 — в Ь(в — 1) в (в2 — в(2Ь2 + 1) + 2Ь2) 0 0
—1
0
0
4
0
0
4(2в — 2Ь2 — 1)
1
2
0
00
00
1
2 0
2 2 в
V — 5 в 5 Ь(в — 1) 10(в2 — в(2Ь2 + 1) + 2Ь2) 0 1 )
01(8,Ь)
—в
—1
0
0
в(Ь — 1) — ^ (в + 2Ь2 — 4Ь + 1)
Ь — 1 — в(2в + 2Ь2 — 4Ь + 1)
0
0
, — в - в(Ь — 1) — — (в + 2Ь2 — 4Ь + 1)
V 5 5 1 ; 101 ;
0
0
0
0
00 Ь — 1 0
00
02(в,Ь)
1 — в Ь(в — 1)
—1
0
0
0
0
в(в2 — в(2Ь2 + 1) + 2Ь2) 0 0
4
4(2в — 2Ь2 — 1)
1
2
0
00
00
Ь0
2 2 в
{ — 5 в 5 Ь(в — 1) 10(в2 — в(2Ь2 + 1) + 2Ь2) 0 1 )
Кроме этого, в нашем случае
0
Ь
Ь
п = 5, г =1, ХТ =( 0 0 0 1 0 ), £15 ^ 0 0 0 1 0 ), 5 = 1,
1 вв(в — 1) + А + Вв ^
[X(в) — Х1 (в)в-1^1^ V+B
— (-в — 1) + В
вв
Т
0
(в — 1) + Вв
+1
4. Исследование нелинейной задачи
Общая методика исследования поставленной задачи основывается на переходе с помощью построенной выше обобщенной матрицы Грина от исходной краевой задачи
23 = А2 + Т(в, 2), М2 = В
(8)
(9)
к некоторой операторной системе и в применении к ней метода последовательных приближений.
Применяя к нелинейной краевой задаче (8), (9) сформулированное выше утверждение и рассматривая нелинейность Т(в, 2(в)) как неоднородность, находим, что задача (8), (9) разрешима тогда и только тогда, когда нелинейность Т (в, 2) £ С [в], С 1[2 ] принадлежит к такому классу вектор-функций, для которых выполнено условие
Р*
0,
(20)
В — М I К(-,Ь)Т(Ь, 2(Ь)) ¿Ь
о
и при этом сводится к однопараметрической системе интегральных уравнений вида
1
2(в, с{) = Х1(в)с1 + [ 0(в, Ь)Т(Ь, 2(Ь, С1)) ¿Ь + [X(в) — Х^в-1в^] V+B. (21)
Заметим, что условие (20) в нашем случае принимает вид
Т3(Ь, 2(Ь, с1)) ¿Ь = 0.
(22)
Рассмотрим численную реализацию решения нелинейной краевой задачи (8), (9) методом последовательных приближений по формуле
1
2(к+1)(в) = Х1(в)с[к+1) + ! 0(в,Ь)Т(Ь,2(к)(Ь) ¿Ь + [X(в) — Х1(в)в-1в15] V+B. (23)
При этом константа с1 определяется единственным образом путем удовлетворения условию (22) на каждом шаге, причем константа с1 добавляется только к и — четвертой компоненте вектора 2. Входящий в формулу интеграл вычисляется по формуле трапеции (с
1
1
о
1
о
о
использованием значений вектора Z с предыдущего итерационного шага). Константа ci из условия разрешимости (22) находится на каждом итерационном шаге при помощи метода Ньютона. В качестве критерия сходимости итерационного процесса используется условие
IIZ(fc+i) - z(k)|| < є.
На начальной стадии (для проверки вида функций, входящих в формулу (23), а также правильного подсчета интеграла) программа была протестирована на линейной неоднородной задаче вида
Zs = AZ + f (s),
MZ = B,
где A, M, B остались без изменения, а вектор правых частей выглядит следующим образом:
( 0 \
1
cos ns
1
1
f (s)
Очевидно, для такого f (в) условие разрешимости (21) выполняется. Задача имеет анали-
тическое решение:
Z
в sin ns Bs(s — 1^ . s(s — 1)
-—З— + ---1 + A + Bs + K J
2
2
в cos ns в(2s — 1^ „ 2s — 1
-—о— + ^—;--------- + B +
2
sin ns
2
+1
V 5
2 в sin ns 2 -------------------+
5
2s — 1 2
РФ — 1) + bs + s(s — 1)
2
2
+ s + 1
/
Поскольку тестовая задача не является нелинейной, метод последовательных приближений к ней не применялся. По результатам теста численное решение совпадает с аналитическим с точностью до 0(10-12).
При решении исходной задачи в качестве начального приближения использовался вектор
Z(0)
( A + Bs\ B
P(s)
C
p(s)
1
(С = 0.001 — 0.004). Функция р имеет вид (см. [5])
р(8)
1, 0 < в ^ — А,
6
11
р1 --А < в < а + А,
6
1
5
5, - + А < в < — — А
6
6
55 р2(в), - — А < в < - + А,
6 - - 6 ’
1,
5
+ А < в < 1,
где
р1(в) р2(в) =
1 — 5
1\3 31 — 5( 1\ 1 + 5
4А3 \в — 6 — Г — 6 +~Г•
1 — 5 (в 5\3 31 — 5 { 5\ 1 + 5
4А3 Г — в) +4~А Г — 6 / + 2 ,
500
А
1
24'
Времена релаксации, входящие в Т, вычисляются по формулам, приведенным в [5]:
т^т = 0.2203 ■ 4 ■ 10-12т
Р,Ы1 1
т = а + Ь (Е — 1) + с е-(1р,™(Е-1)
1 р,ю и,р,'Ш \ ир,'Ш\^ / ' :
6
2
5
где
ар = 0.1153, bp = -0.0068, cp = 0.4988, dp = 1.5137,
aw = 0.4076, bw = 0.0075, cw = 3.1546, dw = 1.4833.
Графически решения при A = 1, B = 1, в = 5 представлены на рисунке.
Как показывают численные эксперименты, на сходимость итерационного процесса
существенно влияют величины в, B и нижняя граница р. При в не выше 102, B < 5 и
р > 0.01 итерационный процесс сходится быстро (около десяти итераций). При изменении параметров в указанных пределах качественный характер поведения решения сохраняется.
Следует отметить, что рассмотренные численные эксперименты носят предварительный характер и проведены с целью проверки правильности работы алгоритма. В дальнейшем планируется более углубленное исследование задачи (8), (9) с применением описанного алгоритма.
Список литературы
[1] GARDNER G. L., Jerome J. W., Rose D. J. Numerical methods for the hydrodynamic device model: subsonic flow. IEEE Trans. Comput.-aided Design, 8, No. 5, 1989, 501-507.
[2] Gardner G. L. Numerical simulation of a steady-state electron shock wave in a submicrometer semiconductor device. IEEE Trans. Electron Device, 38, No. 2, 1991, 392398.
[3] Anile A. M., Maccora C., Pidatella R. M. Simulation of n+ — n — n+ devices by hydrodynamic model: subsonic and supersonic flows. COMPEL, 14, No. 1, 1995, 1-18.
[4] Blokhin a. M., Iohrdanidy a. a., Krymskikh D. a. Numerical investigation of the hydrodynamic model equations for charge transport in semiconductors. Prepr. Sobolev Inst. of Math. Siberian Branch of the Russian Acad. of Sci., №26, Novosibirsk, 1995.
[5] Blokhin A. M., Iohrdanidy A.A., Merazhov I. Z. Numerical investigation of the gas dynamic model equations for charge transport in semiconductors. Prepr. Sobolev Inst. of Math. Siberian Branch of the Russian Acad. of Sci., №33, Novosibirsk, 1996.
[6] Блохин А. М., Крымских Д. А. Численное исследование одной гидродинамической модели переноса носителей заряда в полупроводниках. Математ. моделирование, 9, №3, 1997, 40-50.
[7] Бойчук А. А. Конструктивные методы анализа краевых задач. Наук. думка, Киев, 1990.
Поступила в редакцию 18 ноября 1997 г., в переработанном виде 4 марта 1998 г.
