Интернет-журнал «Науковедение» ISSN 2223-5167 http ://naukovedenie.ru/
Том 8, №1 (2016) http ://naukovedenie. ru/index.php?p=vol8-1
URL статьи: http://naukovedenie.ru/PDF/15TVN116.pdf
DOI: 10.15862/15TVN116 (http://dx.doi.org/10.15862/15TVN116)
Статья опубликована 12.03.2016.
Ссылка для цитирования этой статьи:
Шевченко А.С. Программный комплекс численного решения задачи о баллистическом диоде // Интернет-журнал «НАУКОВЕДЕНИЕ» Том 8, №1 (2016) http://naukovedenie.ru/PDF/15TVN116.pdf (доступ свободный). Загл. с экрана. Яз. рус., англ. DOI: 10.15862/15TVN116
УДК 519.688
Шевченко Алеся Сергеевна
ФГБОУ ВПО «Алтайский государственный университет» Рубцовский институт (филиал), Россия, Рубцовск1 Кандидат физико-математических наук, доцент E-mail: Ibragimova.a.s@mail.ru РИНЦ: http://elibrary.ru/author profile.asp?id=603476
Программный комплекс численного решения задачи
о баллистическом диоде
Аннотация. К настоящему моменту существует довольно много математических моделей, описывающих с той или иной степенью достоверности физические явления в полупроводниковых приборах. В процессе поиска приближенных решений задач физики полупроводников возникает проблема построения численных алгоритмов для этих моделей. Актуальность конструирования подобных алгоритмов не вызывает сомнений, поскольку на сегодняшний день полупроводниковые устройства являются неотъемлемой частью многих электронных приборов.
В данной статье рассматривается, хорошо известная из физики полупроводников, задача о баллистическом диоде и предлагаются некоторые вычислительные алгоритмы, построенные в отличие от обычных конечно-разностных схем, на других идеях. Предложенные алгоритмы основаны на методе установления, применении сглаживающих регуляризующих операторов и идей схем без насыщения.
В качестве основной математической модели взята недавно предложенная гидродинамическая модель переноса заряда в полупроводниках в одномерном случае.
С этой целью написан программный комплекс на языке Object Pascal в среде Delphi 7, реализующий вычислительные алгоритмы.
Приведены подробное описание программного комплекса и результаты численных экспериментов. Полученные результаты являются физически правдоподобными и в целом согласуются с результатами других авторов.
Ключевые слова: программный комплекс; гидродинамическая модель; баллистический диод; стационарные решения задачи; уравнение Пуассона; метод установления; нестационарная регуляризация; схема предиктор-корректор; техника сплайн-функций; интегральные уравнения
1 658225, г. Рубцовск, пр-т Ленина 200-Б 1
Введение. Математическое моделирование физических процессов в полупроводниковых устройствах имеет большое значение в условиях современного стремительного развития микроэлектроники, сопровождающегося возрастающим спросом на повышение функциональности, производительности и улучшение характеристик микросхем и в последнее время превратилось в быстро развивающуюся область прикладной математики.
К настоящему моменту существует достаточно много математических моделей, описывающих физические явления в полупроводниковых приборах. Встает вопрос о конструировании численных алгоритмов и их обосновании для нахождения приближенных решений таких моделей. В качестве математического описания процесса переноса заряда в полупроводниках используется гидродинамическая модель (Anile и Romano [1, 2]), которая с математической точки зрения представляет собой квазилинейную систему уравнений, записанных в форме законов сохранения, и которая получается из системы моментных соотношений для уравнения переноса Больцмана с помощью принципа максимума энтропии.
Постановка задачи. Рассмотрим модель переноса заряда в одномерном случае в безразмерном виде (процесс обезразмеривания подробно описан в работах [3]):
R + J =0,
Jt 2 RE J =RQ+cu J + c121,
(RE )t+ Ix =JQ+cP,
it10 re2 ] =5 REQ+C21J+C221,
(1)
Р = Я -р- (2)
р - электрический потенциал, Я - электронная плотность, Е - энергия электрона,
2
J = Яи, I = Яд - потоки, и - электронная скорость, д - поток энергии, а = — Е -1, Р = Яа,
Q = рх, £> 0 - диэлектрическая постоянная, р = р(х) - плотность легирования (заданная функция на отрезке [0,1]), коэффициенты с, сп, с12, с21, с22 - гладкие функции от энергии Е, выражения для которых приведены, например, в [3, 4].
Для системы (1), (2) сформулируем известную в физики полупроводников задачу о баллистическом диоде. Это одномерная задача, описывающая полупроводниковый прибор, разделенный на три части. Первая и третья области представлены высоколегированным материалом, поэтому они называются п+ областями. В средней же части (в канале) - зона низкой концентрации легировании (так называемая п -область). В качестве характерных величин Ь и N+ будем рассматривать соответственно ширину п+ - п - п+ - канала и плотность легирования п+ - зоны.
Р =
1, если x принадлежи n+ - зоне,
е N+
о =-, если x принадлежи n - зоне,
N
N плотность легирования п - области. Типичный сглаженный профиль функции р(х) изображен на Рис. 1.
<
Рис. 1. Схематическое представление n+ - n - n+ баллистического диода
Граничные условия для системы (1) при х = 0, 1, t > 0, соответствующие известной физики полупроводников задаче о баллистическом диоде следующие:
R(t,0) = R(t, 1) = 1, '
3
Е(*,0) = Е (Г ,1) = -
Начальные данные при I > 0, 0 < х < 1:
Я(0, х) = Я0 (х), ~ 3 (0, х) = 30 (х), Е(0, х) = Е0 (х), /(0,х) = 10(х), ^
Я0(х), Е0(х) > 0.
Краевые условия для уравнения Пуассона (2) при х = 0,1 (? > 0) следующие:
(р{г ,0) = 0, р(М) = В,
где В - напряжение смещения.
(3)
(4)
(5)
При попытке поиска численных решений этой модели мы столкнулись с рядом существенных сложностей, обусловленных нелинейностью уравнений модели и наличием в ней малых параметров. С помощью традиционной теории разностных схем эту задачу решить не удалось.
Задача (1) - (5) в стационарном случае. В стационарном случае математическая модель (1), (2) сводится к смешанной задаче, состоящей из трех уравнений Пуассона, с краевыми условиями на границах х = 0, 1 [5-7]:
г2
йх
й V
йх й 2Х
йх
= Т{а) (2, X, £>, а) = ах22 + а22Х + а32£> + + а502 + Ьса, = (2, X, д,а,Х,р) = -X2 + ^2 2 + Ъ22Х + 632£> + +
(6)
+ Ъ50 + ^^ ~Р) + ПСа,
% = а = 0 при х = 0,1, < ф = 0 при х = 0, ф = В при х = 1.
Компоненты векторов и, q находятся из следующих соотношений:
Р (Е) {0-(1 + а) X - )2},
(7)
и =
q = О(Е) {-0 + (1 + а) X + Оо (Е )2}.
(8)
а
Здесь р = 1, х = 1п Я, 2 = —, X = ^, а1 = -аР(Е)Р0(Е) + Ъ'О(Е)О0(Е) , а йх йх
= -1 + (1 + а) {Ъ'О(Е) - аР (Е)}, а3 = а'¥ (Е) - Ъ О(Е) - ЪР (Е)Р0 (Е), а4 = -Ъ (1 + а) Р (Е),
Ъ2 =(1 + а) {п О(Е) - ш Р (Е)},
а5 = ЪР (Е) Ъз=- 1
(1 + а)
Ъ1 = -ш Р (Е)Р0 (Е) + п О(Е)О0 (Е) 2 , ш Р(Е) - пО(Е) - пР(Е)Р0 (Е) ,
Ъ4 =
1 + а
+ ЪР (Е) Р0( Е),
2 С
а = а( Е) =
5 1 + а 5
- С
11,
Ъ = Ъ( Е) = - с12, ш = ш( Е) = С11 а
5 1 + а
с21 - 5 ес11
1 + а
п = п(Е) =
Ъ = пР (Е) с12 - ъ
22 12 21 ^ 11 ^ Р (Е) =--3-, О( Е) =--3-, Р,(Е) = 1--3
ЕС
5
12
йе^
йе1
, О0(Е) = 1 -
1 + а ес11
с22 ^ ес12
с21 2 ес11
, йа 3 йа йе1 = СС, -СС?, а =-=--и т.д.
11 22 21 12 йа 2 йЕ
Также для исходной задачи (1) - (5) в стационарном случае удалось получить систему интегральных уравнений для определения ф, Е, 1, 3 при Я = р, которую можно воспринимать как приближенную математическую модель для поиска стационарных решений задачи о баллистическом диоде при малых значениях величины а.
Для поиска стационарных решений задачи (1) - (5) с помощью системы интегральных уравнений необходимо:
1. Задать начальные данные для ф, Е . Начальные данные для а рассчитываются по формуле а = 2 Е -1.
1
3
а„( х) = |
2. Вычислить
х
си + сир
Р
-йт
а
коэффициенты: с, сп, с12, с
X X
)(х) = | , а21(х) = |с" + с22^
х
К (х) = | рсайт,
йт.
Р
Р(1 + а) '
а-
х
г( х) = |
22
Р(1 + а)
йт,
й (х) = [«1(1 + а)Р--^ к\ йт,
О I Р Р I
^ = Оц (1) а22 (1)- а_2 (1) а21 (1).
3. Вычислить величины:
й2(х) = И5 (1 + а)Р---к\йт
[ I 2( )Р р(1 + а) |
3 = ■
а.
22
(1)( й, (1)-В)-О.2 (1)1 йг (1)-5 В
IО =
-Оп
(1)(й1 (1)-В) + ап (1)1 йг (1)--В
4. Вычислить
2
2
значения
2
р, а : ю(х) = й(х)-оп(х)3-о12(х)I
3 (а( х ) +1)
а(х) = ®(х)-2й2(х) + 2о21 (х)3 + 2о22 (х)/0, р(х) = ю(х) + а(х), Е(х) =
Пункты 2-4 повторяем до тех пор, пока значения неизвестных не будет удовлетворять заданной точности.
Нестационарная регуляризация краевой задачи (6) - (7). В работах [5-7, 9] для нахождения решения краевой задачи (6-7) использовался метод установления [8]. С этой целью, вместо краевой задачи рассматривалась ее нетривиальная нестационарная регуляризацию Соболева [10], т.е.
(1 =
(9)
т = —, т = — - дифференциальные операторы, переменная ^ играет роль времени.
дХ дх
Замечание. Также рассматривалась нами и параболическая регуляризация [10] уравнений (6):
тр
Ч = %2Х~ Р(ср) (€<7, йх, %<Р,X, р),
но в ходе численных экспериментов было установлено, что скорость сходимости метода установления при использовании регуляризации Соболева (9) намного выше скорости сходимости при применении параболической регуляризации, поэтому в дальнейшем при разработке алгоритмов используем только систему уравнений (9).
Решение системы уравнений (9) ищем в области I > О, О < х < 1 с краевыми условиями (7) и начальными данными
О
О
О
О
О
X
х
ф( а х ) = ф0 ( х ),
2
а(0,х) = а0 (х) = ^Е0 (х)-1, (10)
х( 0, х) = Х> (х) =1п я0 (х ).
Для нахождения решения смешанной задачи (7), (9) - (10) разработано несколько алгоритмов, основанные на сведение задачи к интегральным уравнениям, использование техники сплайн - функций и схемы предиктор - корректор [5-7, 9]. Вкратце опишем эти алгоритмы.
Модельные задачи. Нетрудно заметить, что в системе (9) уравнения различаются только правой частью. Поэтому для построения алгоритмов мы использовали модельную задачу
(1 -£2 )ту = £2 у - /(у), г > 0, 0 < х < 1;
у = 0 при х = 0,1, г > 0; (11)
У = У0 (х) при г = 0, 0 < х < 1,
- известная правая часть. Подставляя в (11) вместо и ф,
а, х вместо у соответственно, получим систему уравнений (9).
Нам потребуется также другая форма записи модельной задачи (11):
\т0+0=/(у), Г>0, 0<х<1;
\ у у ' ' '
У У
< х < .
(12)
= 0у0 (х) =У0 (х) - У0 (х) при г = 0 0 < х <1
еу=еу(их)=у(их)-?у(их), 7(>,) = /(>,)(г,х) =у(г,х)-/(>,)(г,х).
Для нахождения стационарных решений (1) - (5) с помощью метода установления в уравнении (1 -£2 )ту = £2у - /(у) проведём дискретизацию по переменной Для этого введём
обозначения: уп(х) = у(пА,х) = у, п = 0,1,..., У = уп+1 (х), А - шаг разностной сетки по времени Аппроксимируя производную ту выражением У-у получаем:
£2У = ВУ + в = —, £ = ^ у у + А^(у). (13)
1+ А 1+ А
Техника сплайн-функций. Приближённое решение модельной задачи (11) ищем в виде интерполяционного кубического сплайна класса С2 (см. [11]):
и1
S(х) = (1 -т)Ук + ТУк+1 - — т(1 -т)[(2-т)ш^ + (1 + т)шк+ ], (14)
6
х х<
т=—^, х е[хк, Хк+1 ], к = 0, N -1, хк = к • к , N • И = 1, у = У(хк), шк = ^У(хк).
к
Первая производную от кубического сплайна следующая:
£У(х) = - к [(2 - 6т + 3г2М + (1 - 3т2)тк+1 ], х е^, х^+1 ], к = О, N -1.
Вычисляя ^(х* + О), ^(х* - О), где ^(хк + О) = 7+1 - 7 - — [2т, + тк+1 ].
к 6
У - 7 к
^^(х, - О) = —-—--[т, ч + 2щ ] и приравнивая их, получаем
к 6
^1
1 1 3 _ -
"тк-1 + 2тк +-тк+1 = 77 (Лук-Лук ), к =1N -1, (15)
2 2 к
/ = щ-1, // = 1 - - разностные операторы, , операторы сдвига:
7к = Ук±1 1 = .
Полагая в (13) х = х, и подставляя ; Ук в (15), получаем следующую трехслойную
разностную схему:
= {&-1 + + &+1} ,к = 1, N -1, & = & (хк).
Краевые условия при к = О и к = N для нахождения функции Ук ,к = О, N с помощью системы уравнений (16) следуют из краевых условий для уравнения (11):
уо = О,УN = О. (17)
Систему (16) - (17) решаем методом прогонки [8].
Замечание. Значения производной ;у, которые потребуются для расчётов, в узлах разностной сетки (т.е. при х = кк = хк, к = О, N ), вычисляются с помощью сплайна (14).
Так как ^(х) = - к [(2 - 6т+ 3г2К + (1 - Зг2^ ], х е^, хк+х ], к = О, N -1,
то ^(хк+1) = ^^-^[тк + 2тк+1 ], х е [хк,хк+1 ], к = О^, £У(0) = ^ + ^[2то + т]. к 6 к 6
Сведение модельной задачи (12) к системе интегральных уравнений. Задача (12) -неоднородное дифференциальное уравнение (с параметром х). Решение данного уравнения следующее:
Х
ву (х, х ) = в~вуо (х) +[ у (С, х)-/у)(С, х) Ус, (18)
о
где вуо (х) =у0 (х) -;2у0 (х) - начальные данные для соотношения (12).
Из ву (х,х) =у(х,х)-;2у (х,х) и у = О при х = 0, 1, Х > 0 получаем краевую задачу с параметром X:
(у),
у (г, х)- у (г, х) = -0у (г, х), [у (г,0) = у (г,1) = 0.
Решение полученной краевой задачи (19) запишем так:
(19)
1
у (г, х) = - [ О (х, ^)ву (г, ^) ж,
0
(20)
О ( х, ^) = ^
бЬ ^ бЬ ( х -1)
бь х бь (^ -1) бь1
0 < ^ < х,
функция Грина.
, х < ^ < х.
Дифференцируя (20) по х, получаем выражение для £у, которое будет использовано нами для расчётов правых частей уравнений в системе (9):
х сИх 1
£у (г, х) = -[ Ск (х - 5)0ф(г, 5уъ —— | Ск (5 -1) вф (г, 5й.
о о
(21)
Численные расчеты для модельной задачи (11) проводятся с переменным шагом И = хг+1 - х, i = 0, N-1 на сетке по х. В интегральных уравнениях (18), (20), (21) проводим
дискретизацию по переменной I и вводим обозначения: (в). =в(пА,х), у. = у (пА,х ), А -шаг разностной сетки по переменной В результате переходим к соотношениям:
1
у. =-{ О ( х, 5 ву ((п - 1)А, 5 ) Ж,
£у(пА,х) = -] сИ(х -5)вф((п- 1)А,5й-^}Ск(5- 1)вф((п- 1)А,5(22)
0 0
(ву )"
(
= е
-пА
пА
в( 0, х )+| е^ у (С, х )-х )
Л
йС
ву (0,х ) =у0 (х)-£2у, (х), п = 1,2,....
Интегралы в системе (22) вычисляем методом трапеции.
Схема предиктор-корректор. Нам понадобятся уравнение (12) и интегральные соотношения, полученные выше, в которых проводим дискретизацию по переменной
Для уравнения (12) схема предиктор - корректор запишется так:
в-в +в=
у
И"-
/ЛП+\ /лп
И'-
А/2 ' V / А Выражая в** из первого, а в'П+1 - из второго уравнения, получаем
0
в;+А(/(у)]г
Л* _ } 2У / йп+1_ ву =-Д-, ву ='
1+-
1 + д
2
Используя уравнения (20) - (21), переходим к следующей системе выражений:
(в)"+-(/м (/{у>)" = у?-/{у>{"&,*<). (*;),=— а
у* = -[а (х, 5 )в* (5) ^,
о
кг
1
уГ1 = -[ а (х, 5 ву ((п+1)д, 5 ) ^
ы
(23)
1 + д
;у (( п + 1)д, х- ) = - )ск ( х - 5 )вр(( п + 1)д, 5 у - ^} ск ( 5 - 1)вД( п + 1)д, 5
о о
Интегралы в системе (23) вычисляем методом трапеции.
Модификация алгоритмов. Теперь рассмотрим модификации алгоритмов. С вычислительной точки зрения, перспективным оказалось использование для расчёта потенциала р и его производной следующих соотношений:
рп = р) а ( х-, 5 )(п-1)д* )-Р( 5 ))
;р( пД, х- ) = ^[( е
-Р( 5 )У-^}(1 - 5 )(е
-р(
5 ) Ш
а (х, 5 ) =
|5 (х -1) ,0 < 5 < х, I х (5 -1), х < 5 < х.
Решение регуляризованной системы с помощью описанных методов.
Для осуществления расчётов нужно определить начальные значения переменных (Е, Я , р ),а затем выполнить следующие шаги:
1. В начале вычисляем Т{ср)(х'\ Р), где х"=х{п^,х) ~ значение X с предыдущего слоя. Затем находим решение ри+1 =р(( п + 1)Д, х) первого уравнения задачи (9) и его производную ;рп+х на следующем слое.
2. Рассчитываем величину Т{а) (%с> (иА, х), & («4 *), й<Р ((<п +1) А, *) + Я, ег" ). Далее решаем второе уравнение задачи (9). Находим аи+1 =а(( п + 1)Д, х) ,;ап+1.
л
2
*
<
о
о
3. Вычисляем 0={х) (¿¡а((и + 1)А,х),&(пА,х),£<р((и +1) А,х) + В,ап+1 ,х",р)-
Далее решаем третье уравнение задачу (9). Находим = х((п + 1)А, х
Выполнив 1-3 шаги, алгоритм переходит на следующий временной слой. Данные шаги выполняются до тех пор, пока не получим решение заданной точности.
При реализации вычислительных алгоритмов, мы использовали идею метода установления и проводили итерации по временной переменной. В ходе вычислений возникал скачкообразный рост при некоторых значений переменных задач, обусловленный ее нелинейностью: норма решений для энергии становилась большой, что вызывало переполнение буфера и остановку работы программы. Для решения этой проблемы применили итерации по нелинейности.
Основная идея алгоритма, базирующегося на итерациях по нелинейности, состоит в том, чтобы пересчитать значения параметров и переменных задач по формулам, предназначенным для перехода на следующий временной слой, но оставаться при этом на текущем слое.
Опишем, например, технику сплайн-функций, с применением техники итераций по нелинейности.
Трехслойная схема (16) преобразуется к виду:
1 - h- ©| й- 2 ji+h2 ©I Yw+ji - ©I Y;1] =
(24)
12 _ _
h-{ЗЁ-Ч+ 4Fl;-1]+ Fl:r1]}, k = 1,N-1, ; = 1,M.
Здесь М — количество итераций по нелинейности на каждом временном слое, значения 1, к = 1, N -1, вычисленные на (т — 1) - й итерации по нелинейности.
Алгоритм работы схемы следующий. Значения , к = 1, N — 1 берутся с предыдущего
временного слоя. На ; - й итерации по нелинейности вычисляются значения , k = 1, N -1, затем по этим значениям по (24) вычисляются Fk"^ и программа переходит на (; +1) -у
итерацию по нелинейности. Как только ; = M программа переходит на следующий временной слой.
Описание программного комплекса. Для нахождения приближенных решений в задачи о баллистическом диоде, написан программный комплекс на языке Object Pascal в среде Delphi 7 [12], реализующий описанные выше вычислительные алгоритмы.
Программный комплекс состоит из следующих модулей:
Unit_main.pas - основной модуль программы.
Unit_matrix.pas - содержит процедуры и функции: сложение, вычитание и умножение двух матриц; умножение матрицы на число; умножение матрицы на вектор; сложение двух векторов; вычисление определителя и обратной матрицы; методы Гаусса и прогонки для решение систем линейных алгебраических уравнений; вычисление первых производных с помощью кубического сплайна класса С2, вычисление вторых производных.
Unit_c.pas - содержит значения физических параметров для баллистического диода; константы электронно-фононного взаимодействия, энергии фононов; функции для
вычисления обычных и модификационных функций Бесселя; процедуры для вычисления коэффициентов с,с , с ,с , с и их производных; процедуры для вычисления параметров, представленных выше алгоритмов.
Unit_decision.pas - содержит процедуры реализации предложенных алгоритмов.
Главная форма программного комплекса содержит меню из 3 пунктов: «Файл», «Численное решение», «Справка» (см. рис. 2).
Меню «Файл» содержит подменю: «Считать данные из файла», «Сохранить в файл», «Выход».
.Ф Численный анализ задачи о баллистическом диоде
Файл Численное решение Справка
Считаггь данные из файла
Сохранить в файл ► Результаты расчета Графики
Вьиод
Рис. 2. Главная форма программного комплекса
При выборе пункта меню «Численное решение» появляется окно «Численное решение задачи о баллистическом диоде», см. рис. 3.
Рис. 3. Окно расчета численного решения задачи о баллистическом диоде
http://naukovedenie.ru 15туШ16
Чтобы найти стационарные решения задачи о баллистическом диоде необходимо в программе выполнить несколько шагов.
Шаг 1. Задать параметры:
а) физические параметры, характеризующие исходную задачу;
б) численные параметры, определяющие вычислительный процесс;
в) начальные данные для искомых функций р , Е, Я , 1, 3 .
Начальные значения переменных задачи лучше задавать равным значениям этих переменных в состоянии глобального термодинамического равновесия.
Необходимые для расчетов физические параметры и их описание содержаться в табл.
1.
Таблица 1
Значение физических параметров
Величина Описание Значение
N+ Плотность легирования п+ - области 5 х 1017 \ ст
N Плотность легирования п - области 2 х1015 \ ст
Ь Длина п+ - п - п+ - канала 6 х10-7 т или 3х10-7 т
Р 10749.32 или 26873.33
8 Значение функции р( х) в п -зоне баллистического диода ^ т 1.19170 х105 —
V Напряжение смещения 0.5 ^ 2Volt
Функция р в расчетах является кусочно-линейной версией функции плотности легирования р( х) (см. рис. 1) и задается в программе так:
р( х н
15
1, если 0 < х < ё , Р ( х), если ё < х < ё + ^, 8, если 1 - ё < х < 1 - ё + ^, 1, если 1 - ё + 4 < х < 1,
Р (х) = -Тв -8)
1 2 „з „
— X" —X + X
5 3
1 + 3 _ х-ё
л--,х =-
2 А,
Р (х ) = 15 (1 -8)
1 2 .з _
— х~ —X + X
5 3
\ + 8 „_х-(\-ё)
Н , X — .
2 А,
Безразмерная величина В вычисляется автоматически, как только задали значение напряжение смещения V.
В табл. 2 приведены необходимые для расчетов численные параметры и диапазоны их значений.
Алгоритмы прекращают работу, если достигается необходимая точность £1:
N /
£( - Е?
г=0
_п+1 _п
а1 -а1
+
п+1 п
Р -Р
< &
Таблица 2
Значение численных параметров
Величина Описание Значение
N Количество узлов сетки по х 100 ^ 500
М Шаг разностной сетки по 1 0.02 -И0-5
м Количество итераций по нелинейности 1^10
й 1 1 — или — 6 3
& Точность вычислений
4 Напряжение смещения 1 1 8 ' 12
После того как указали все необходимые параметры необходимо нажать кнопку «Утвердить начальные данные».
Шаг 2. Выбор алгоритма решения задачи (интегральные уравнения, модификация интегральных уравнений, схемы предиктор-корректор, техника сплайн-функций, приближенная математическая модель), т.е. нужно выбрать соответствующую вкладку с алгоритмом решения задачи.
Управление расчетом осуществляется с помощью кнопок «Начать расчет» и «Прервать расчет» на выбранной вкладке.
Кнопка «Начать расчет» становится доступной только после того как будет нажата кнопка «Утвердить начальные данные».
При нажатии кнопки «Начать расчет», ее название меняется на «Идет расчет».
Кнопка «Прервать расчет» становится доступной только после того как будет нажата кнопка «Начать расчет».
Все полученные результаты отражаются на вкладках: «Данные расчета», «График плотности», «График зависимости энергии Е от Ь>, «График приближенного решения для потенциала», «График приближенного решения для энергии Е», «График приближенного решения для электронной плотности Я», «График приближенного решения для I», «График приближенного решения для У», «График приближенного решения для и» и «График скорректированного и». Только для алгоритма «Приближенная математическая модель» есть вкладка «Вольт-амперная характеристика баллистического диода».
Вкладка «График зависимости энергии Е от Ь> (наибольшее значение энергии на п-ом временном слое) позволяет нам следить за поведением энергии (см. Рис. 4).
Кнопки «Сохранить результаты расчета в файл» (пункт меню «Файле» ^ «Сохранить в файл» ^ «Результаты расчета»), «Считать данные из файла» (пункт меню «Файле» ^ «Считать данные из файла») позволяют осуществлять сохранение и загрузку результатов расчета (приближенные решения для электрического потенциала, плотности электронов, энергии электронов, электрического тока и потока энергии). Данные сохраняются в файлах с расширением .Ш, что позволяет использовать полученные результаты в качестве начальных условий для дальнейших расчетов.
Графики полученных решений можно сохранить в любом графическом формате с помощью панели «Сохранение графиков» или пункта меню «Файле» ^ «Сохранить в файл» ^ «Графики».
Рис. 4. Вкладка «График зависимости энергии Е от I»
Результаты тестовых расчетов. Приведем теперь графики численных решений, полученные в результате расчетов.
На рис. 5. приведены графики энергии, потенциала и вольт-амперная характеристика баллистического диода (зависимости силы тока от напряжения смещения В).
Рис. 5. Результаты расчетов для случая £ = 0 с параметрами V = \5Volt, а = 0.004,
Ь = 6 х10-7 т, N = 100, £ = 10-6
Расчеты проводились с помощью различных алгоритмов.
Вычислительный алгоритм для системы интегральных уравнений и его модификация дают хорошие результаты при малом а = 0.004 _
Приближенные решения, полученные с помощью модификации интегральных уравнений, представлены на Рис. 6.
0,000 0,113 0,225 0,325 0,433 0,542 0,650 0,762 0,863 0,975 0,000 0,113 0,225 0,325 0,433 0,542 0,650 0,762 0,863 0,э75
x .v
Рис. 6. V = 1.5 Volt, < = 0.004, L = 3х10"7m, N = 100, Д1 = —, е1 = 10"6
1 12 1
Другие подходы (с использованием техники сплайн-функций и схемы предиктор-корректор) дают хорошие результаты при не больших <.
В результате численных экспериментов было выявлено, что скорость сходимости алгоритма, использующего технику сплайн-функций выше, чем у алгоритма, использующий схему предиктор-корректор.
C целью верификации предложенных алгоритмов для задачи (9) проводились расчеты с использованием метода ортогональной прогонки, который позволяет в принципе проводить расчеты с гарантированной точностью. На рис. 7. приведены результаты сразу трех алгоритмов (схема предиктор-корректор, техника сплайн функций и метод ортогональной прогонки) и видно, что результаты практически совпадают.
-- 1 ----II
0,000 0,080 0,160 0,240 0,520 0,400 0,480 0,560 0,640 0,720 0,800 0,880 0,960 0,000 0,080 0,160 0,240 0,320 0,400 0,480 0,560 0,640 0,720 0,800 0,880 0,960
* х
Схема предиктор-корректор Метод ортогональной прогонки 'Техника сплайн-функций
Рис. 7. V = \Volt, а = 0.2, Ь = 6х10-7т, N = 200, А, = —, £ = 10-6
1 12 1
Заключение. Проведено множество численных расчетов. Полученные результаты являются физически правдоподобными и в целом согласуются с результатами других авторов. Данные результаты носят как теоретический, так и вычислительный характер. Данная работа представляет интерес для специалистов в области вычислительной математики, может быть использована для моделирования физических процессов в полупроводниковых устройствах и решения различных прикладных задач.
ЛИТЕРАТУРА
1. Anile A.M., Romano V. Non parabolic band transport in semiconductors: closure of the moment equations // Cont. Mech. Thermodyn, 1999. - V.11. - P. 307-325.
2. Romano V. 2D simulation of a silicon MESFET with a non-parabolic hydrodynamical model based on the maximum entropy principle // J. Comp. Phys, 2002. - V.176. - P. 70-92.
3. Blokhin A.M., Bushmanov R.S., Romano V. Nonlinear asymptotic stability of the equilibrium state for the MEP model of charge transport in semiconductors, // Nonlinear Analysis, 2006. - V. 65. - P. 2169-2191.
4. Blokhin A.M., Ibragimova A.S. Numerical method for 2D Simulation of a Silicon MESFET with a Hydrodynamical Model Based on the Maximum Entropy Principle // SIAM Journal on Scientific Computing, - 2009. - V.31. - Issue 3. - P. 2015-2046.
5. Blokhin A.M., Ibragimova A.S. 1D Numerical Simulation of the MEP Mathematical Model in ballistic diode problem // Journal of Kinetic and Related Models, - 2009. -V. 2. - N.1. - P. 81-107.
6. Блохин А.М., Ибрагимова А.С., Семисалов Б.В. Конструирование вычислительных алгоритмов для задачи о баллистическом диоде // Вычислительная математика и математическая физика, - 2010. - Т.50. - №1. - c. 1-21.
7. Blokhin A.M., Ibragimova A.S., Semisalov B.V. Design of Numerical Algorithms for the Ballistic Diode Problem // Computational Mathematics and Mathematical Physics, - 2010. - V. 50. - N. 1. - P. 180 - 200.
8. Бабенко К.И. Основы численного анализа. - Москва-Ижевск: НИЦ "Регулярная и хаотическая динамика", 2002. - 848 c.
9. Блохин А.М., Ибрагимова А.С., Семисалов Б.В. Численный анализ задач переноса заряда в полупроводниковых устройствах. Saarbrucken, Germany: Palmarium Academic Publishing, 2012. - 209 c.
10. Шевченко А.С. Алгоритм поиска приближенных решений уравнения Пуассона // Молодой ученый, 2015. - №3. - С. 18-23.
11. Завьялов Ю.С., Квасов Б.И., Мирошниченко В.Л. Методы сплайн-функций. -М.: Наука, 1980. - 355 c.
12. Белов В.В., Чистяков В.И. Программирование в Delphi: процедурное, объектно-ориентированное, визуальное. Учебное пособие для вузов. - М.: Горячая линия-Телеком, 2014. - 240 с.
Shevchenko Alesya Sergeevna
Altay State University Rubtsovsk branch, Russia, Rubtsovsk E-mail: Ibragimova.a.s@mail.ru
Software complex of numerical solution of the ballistic
diode problem
Abstract. At present there are quite a few of mathematical models describing physical phenomena in semiconductor devices with one or another degree of reliability. The problem of constructing numerical algorithms for these models arises in process of finding approximate solutions of semiconductor physics. Undoubtedly, the designing of such algorithms is actual, because today semiconductor devices are essential parts of many electronic devices.
In this paper, we consider the ballistic diode problem, which is well known in the physics of semiconductors, and propose computational algorithms based on ideas other than those underlying usual finite-difference schemes. The proposed algorithm is based on the stabilization method, the application of regularized smoothing operators and ideas of schemes without saturation.
As the basic mathematical model, we use a recently proposed hydrodynamic model of charge transport in semiconductors in one-dimensional case.
The computational algorithms are implemented as a software package using Object Pascal in Delphi 7 environment.
The detailed description of a software complex and results of numerical experiments are provided. The obtained results are physically plausible and, in general, agree with the results of other authors.
Keywords: software complex; hydrodynamic model; ballistic diode; stationary solutions of problem; the Poisson equation; stabilization method; nonstationary regularization; predictor-corrector scheme; spline functions technique; integral equations
REFERENCES
1. Anile A.M., Romano V. Non parabolic band transport in semiconductors: closure of the moment equations // Cont. Mech. Thermodyn, 1999. - V.11. - P. 307-325.
2. Romano V. 2D simulation of a silicon MESFET with a non-parabolic hydrodynamical model based on the maximum entropy principle // J. Comp. Phys, 2002. - V.176. - P. 70-92.
3. Blokhin A.M., Bushmanov R.S., Romano V. Nonlinear asymptotic stability of the equilibrium state for the MEP model of charge transport in semiconductors, // Nonlinear Analysis, 2006. - V. 65. - P. 2169-2191.
4. Blokhin A.M., Ibragimova A.S. Numerical method for 2D Simulation of a Silicon MESFET with a Hydrodynamical Model Based on the Maximum Entropy Principle // SIAM Journal on Scientific Computing, - 2009. - V.31. - Issue 3. - P. 2015-2046.
5. Blokhin A.M., Ibragimova A.S. 1D Numerical Simulation of the MEP Mathematical Model in ballistic diode problem // Journal of Kinetic and Related Models, - 2009. -V. 2. - N.1. - P. 81-107.
6. Blokhin A.M., Ibragimova A.S., Semisalov B.V. Konstruirovanie vychislitel'nykh algoritmov dlya zadachi o ballisticheskom diode // Vychislitel'naya matematika i matematicheskaya fizika, - 2010. - T.50. - №1. - c. 1-21.
7. Blokhin A.M., Ibragimova A.S., Semisalov B.V. Design of Numerical Algorithms for the Ballistic Diode Problem // Computational Mathematics and Mathematical Physics, - 2010. - V. 50. - N. 1. - P. 180 - 200.
8. Babenko K.I. Osnovy chislennogo analiza. - Moskva-Izhevsk: NITs "Regulyarnaya i khaoticheskaya dinamika", 2002. - 848 c.
9. Blokhin A.M., Ibragimova A.S., Semisalov B.V. Chislennyy analiz zadach perenosa zaryada v poluprovodnikovykh ustroystvakh. Saarbrucken, Germany: Palmarium Academic Publishing, 2012. - 209 c.
10. Shevchenko A.S. Algoritm poiska priblizhennykh resheniy uravneniya Puassona // Molodoy uchenyy, 2015. - №3. - S. 18-23.
11. Zav'yalov Yu.S., Kvasov B.I., Miroshnichenko V.L. Metody splayn-funktsiy. -M.: Nauka, 1980. - 355 c.
12. Belov V.V., Chistyakov V.I. Programmirovanie v Delphi: protsedurnoe, ob"ektno-orientirovannoe, vizual'noe. Uchebnoe posobie dlya vuzov. - M.: Goryachaya liniya-Telekom, 2014. - 240 s.