автореферат диссертации по информатике, вычислительной технике и управлению, 05.13.18, диссертация на тему:Математическое моделирование в задачах переноса заряда в полупроводниковых кремниевых устройствах

На правах рукописи

СЕМИСАЛОВ Борис Владимирович

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

05.13.18 - математическое моделирование, численные методы и комплексы программ

АВТОРЕФЕРАТ

диссертации на соискание ученой степени кандидата физико-математических наук

Омск - 2011

2 6 МАЙ 2011

4848019

Работа выполнена на кафедре дифференциальных уравнений Новосибирского государственного университета

Научный руководитель: доктор физико-математических наук,

профессор Блохин Александр Михайлович

Официальные оппоненты: доктор физико-математических наук,

профессор Воеводин Анатолий Фёдорович

доктор физико-математических наук, профессор Задорин Александр Иванович

Ведущая организация: Институт вычислительных технологий СО РАН

Защита состоится 16 июня 2011 года в 14:30 на заседании объединенного совета по защите докторских и кандидатских диссертаций ДМ 212.179.07 при Омском государственном университете им. Ф.М. Достоевского по адресу: 644099, г. Омск, ул. Певцова, 13.

С диссертацией можно ознакомиться в библиотеке Омского государственного университета им. Ф.М. Достоевского по адресу: 644077, г. Омск, Пр. Мира, 55-А.

Автореферат разослан « /2 » мая 2011 г.

Ученый секретарь диссертационного совета

к.ф.-м.н., доцент

■ШД-

А.М. Семёнов

ОБЩАЯ ХАРАКТЕРИСТИКА РАБОТЫ

Актуальность темы. На фоне стремительного развития микро- и на-ноэлектроники и уменьшения размеров полупроводниковых устройств появляется множество новых математических моделей, описывающих с той или иной степенью достоверности физические явления, протекающие в этих устройствах. В процессе поиска приближённых решений задач физики полупроводников возникает проблема построения численных алгоритмов для указанных моделей. Актуальность конструирования подобных алгоритмов не вызывает сомнений, поскольку на сегодняшний день полупроводниковые устройства являются неотъемлемой частью многих электронных приборов. От характеристик этих устройств существенно зависит мощность и надёжность современных вычислительных машин. Рассчитать все характеристики полупроводника с необходимой точностью - и есть задача математического моделирования.

Моделирование процесса переноса зарядов в полупроводниках субмикронного размера основывается на кинетическом уравнении Больцма-на для функции распределения носителей зарядов. Однако прямое численное интегрирование этого уравнения с помощью метода Монте-Карло требует существенных вычислительных затрат и в некоторых случаях (например, если концентрация носителей заряда в отдельных областях полупроводникового устройства мала) не позволяет получить точный результат. Разумный компромисс между физической точностью и вычислительной эффективностью может быть достигнут посредством выделения моментов уравнения Вольцмана. Простейшей моделью переноса заряда, полученной подобным методом, является дрейф-диффузионная модель Schockley и van Roosbroeck'a1, которая впоследствии модифицировалась, дополнялась и усложнялась многими исследователями. Она состоит из уравнений неразрывности для концентраций носителей зарядов (электронов и дырок) и уравнения Пуассона для электрического потенциала. Однако уменьшение размеров полупроводниковых устройств требует расширения общепринятой дрейф-диффузионной модели, с целью учёта энергии носителей зарядов. Эта цель достигается в гидродинамических моделях. Одна из первых гидродинамических моделей представлена в работе Blotekjaer'a2 и содержит законы сохранения для числа частиц, импульса и энергии носителей зарядов, а также уравнение Пуассона для электрического потенциала.

В диссертации рассматривается недавно предложенная гидродинамическая МЕР модель3. Эта модель представляет собой квазилинейную систему

^■¡.oosbroek W. van. Theory of flow of electrons and holes in germanium and other semiconductors // Bell System Techn.J. 1950. Vol. 29. P. 560-607.

2Blotekjaer K. Transport equations for electrons in two-valley semiconductors // IEEE Trans. Electron Devices. 1970. Vol. ED-17. P. 38-47.

3См. работу: Anile A. M., Romano V. Non parabolic band transport in semiconductors: closure of the moment equations // Cont. Mech. Thermodyn. 1999. Vol. 11. P. 307-325.

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

Как мы видим, математические модели полупроводниковых устройств эволюционируют и усложняются. Большинство из них содержат нелинейные уравнения в частных производных. Таким образом, помимо проблемы конструирования вычислительных алгоритмов, остро встаёт вопрос об их обосновании. Разрешить этот вопрос можно лишь в ходе детального анализа соотношения между исходной математической моделью явления и её разностным аналогом. В связи с этим, используя принципы монографии4, при конструировании и исследовании вычислительных моделей в диссертации мы отталкивались от требования адекватности вычислительной модели исходной дифференциальной задаче. Под адекватностью понимаем следующее: разностная модель строится так, чтобы с её помощью можно было доказать теорему существования решения исходной задачи. Подобные доказательства опираются на наличие разностных аналогов априорных оценок, которые удалось получить для уравнений МЕР модели в некоторых случаях (при дополнительных предположениях). Последнее обстоятельство представляется чрезвычайно важным фактом, поскольку благодаря наличию указанных оценок при проведении вычислений мы можем быть уверены в том, что приближённое решение действительно удовлетворяет свойствам исходной задачи.

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

Чтобы преодолеть указанные трудности, желательно организовать дискретизацию таким образом, чтобы происходила автоматическая адаптация позиций узлов сетки к свойствам решения, погрешность аппроксимации была минимальной, а алгоритм, её использующий, давал высокую точность при небольшом количестве узлов. Подобные требования легли в основу вычислительного алгоритма, предложенного для поиска решений уравнений МЕР модели в гл. 1 диссертации. Вслед за автором описанной идеи, К.И. Бабенко будем именовать схемы, автоматически учитывающие априорные

4 Блохин A.M., Алаев Р. Д. Интегралы энергии и их приложения к исследованию устойчивости разностных схем. Новосибирск, 1993.

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

Основываясь на идеях схем без насыщения и принципе адекватности в диссертации предложена эффективная и удобная технология построения вычислительных моделей для поиска решений 1D и 2D нестационарных смешанных краевых задач. В Приложениях III и IV данная технология успешно применена при поиске численных решений задач, связанных с аэростроением и добычей нефти, а именно при исследовании вопроса об устойчивости ударных волн в сжимаемом вязком газе и при изучении явления параметрического резонанса в слоистых газосо-держащих структурах. Таким образом, построенная технология открывает широкие возможности для поиска приближённых решений различных актуальных прикладных задач.

Цели работы:

1. Конструирование и обоснование вычислительных моделей для поиска стационарных решений уравнений гидродинамической МЕР модели в 1D и 2D случаях.

2. Создание на базе построенных моделей работоспособных алгоритмов и их реализация на ЭВМ.

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

Объектами исследования являются:

- гидродинамическая МЕР модель, описывающая процесс переноса заряда в полупроводниковых приборах;

- задачи о баллистическом диоде и о переносе заряда в транзисторе MOSFET, поставленные в главах работы для уравнений МЕР модели;

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

Основные результаты, выносимые на защиту:

1. Разработаны и теоретически обоснованы два новых подхода для поиска стационарных численных решений уравнений гидродинамической МЕР модели:

- первый подход является оригинальным сочетанием набора регуляриза-ций, идей схем без насыщения, сплайн-интерполяций и метода установления;

- второй подход представляет новый способ конструирования и исследования устойчивых разностных схем для нелинейных задач.

0 Бабенко К. И. Основы численного анализа. М.; Ижевск: НИЦ «Регулярная и хаотическая динамика», 2002.

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

3. При использовании технологии подхода 1 разработан вычислительный алгоритм для поиска стационарных решений задачи о переносе заряда в транзисторе МОБРЕТ.

4. Все алгоритмы, описанные в работе, реализованы на ЭВМ, проведены численные эксперименты, результаты которых представлены в диссертации.

Научная новизна работы. Методы исследования. Автор видит новизну полученных результатов в следующем.

Во-первых, в диссертации предложено два новых подхода для поиска стационарных численных решений уравнений МЕР модели. Структура и методы указанных подходов отражены на рис. 1.

Регулррлзсал^нэя

'неегацмнарнм систем!

1— Метод прямых Система ОДУ второго порйдна

тт

Интерпезкши сплайнами

Трехточечная разностная схема.

•тод прогонки

Решение регуляриэоеанных уравнения на каждом временном слое

ДИСКРСТЯМПКИ

Класс "устойчивых" разностных схем

ТГ

Трёхточечная матричная схема

Метод матричной прогонки ^ "

Решение на каждом «ременном слое

ж;

Метод установления

НЕ

Стэциойарные решения

Рис. 1. Два подхода для поиска стационарных решений уравнений МЕР модели

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

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

В Ш случае: аппроксимация производной неизвестной функции по времени разностным отношением —> сведение задачи к поиску решения ОДУ второго порядка с краевыми условиями —> использование сплайн-функций для интерполяции решения полученного ОДУ —> сведение задачи к трёхто-чечной схеме с граничными соотношениями —> поиск решения на каждом шаге по времени методом прогонки.

В 2В случае: аппроксимация производной неизвестной функции по времени разностным отношением, а производной по одному из пространственных направлений - интерполяционным многочленом с узлами интерполяции в нулях полинома Чебышева —» сведение проблемы к краевой задаче для системы ОДУ второго порядка —» поиск решения полученной системы в виде интерполяционного кубического сплайна класса С2 —* вывод трёхточечных матричных соотношений с граничными условиями —> поиск решения на каждом временном слое методом матричной прогонки.

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

Теоретическая и практическая значимость.

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

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

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

4. Техника, использованная в диссертации для вывода априорных оценок и их разностных аналогов, представляет определённый интерес и может быть применена при исследовании различных задач и уравнений в частных производных.

Апробация работы. Основные результаты диссертации были представ-

лены на: XLVII, XLVIII Международных научных студенческих конференциях «Студент и научно-технический прогресс» (Новосибирск, 2009, 2010 гг.); в рамках международной школы-семинара «International School and Seminar on Modern Problems of Nanoelectronics, Micro- and Nanosystem Technologies Proceedings "INTERNANO-2009»> (Новосибирск, 2009 г.); международной конференции «International Conference on Computational Technologies In Electrical and Electronics Engineering "Sibircon-2010>> (Иркутск, Листвянка, 2010 г.); на III международном форуме по нанотехнологиям (Москва, 2010 г.). Кроме того, результаты работы докладывались на семинарах: в Ин-те динамики систем и теории управления, г. Иркутск (председатель д-р физ.-мат. наук A.B. Синицын); «Теоретические и вычислительные проблемы задач математической физики» в Ин-те математики СО РАН (рук. проф. A.M. Блохин); «Общеинститутский математический семинар» в Ин-те математики СО РАН в рамках конкурса Сиб. мат. журнала; «Информационно-вычислительные технологии» в ИВТ СО РАН (рук. акад. Ю.И. Шокин, проф. В.М. Ковеня); «Прикладная гидродинамика» в Ин-те гидродинамики СО РАН (рук. чл.-корр. РАН В.В. Пухначёв); «Математическое моделирование и вычислительные методы» в Омском филиале Ин-та' математики СО РАН (рук. проф. А.И. Задорин).

Публикации. По теме диссертации опубликовано 11 работ, куда входят: 4 статьи в изданиях, рекомендованных ВАК, 2 - в зарубежных журналах, 3 - в трудах международных конференций и семинаров, 2 - в тезисах конференций.

Личный вклад автора. При подготовке работы [1] соискатель занимался конструированием и применением алгоритма для поиска стационарных решений уравнений гидродинамической МЕР модели. Автором разработаны и реализованы на ЭВМ все описанные в [2], [4] вычислительные алгоритмы, выполнена дополнительная модернизация алгоритмов с целью повышения их эффективности и обеспечения сходимости метода установления. Работы [3], [5] выполнены при равном вкладе авторов. Соискателем разработаны все вычислительные схемы работы [6] и написан комплекс программ на их основе.

Структура диссертации. Диссертация состоит из введения, четырёх глав, заключения, списка литературы и четырёх приложений. Объём работы - 184 страницы. В диссертации содержатся 43 рисунка и 9 таблиц. Список литературы состоит из 105 источников.

Автор выражает глубокую благодарность и признательность своему научному руководителю д-ру физ.-мат. наук, профессору Александру Михайловичу Блохину за постоянное внимание и помощь в работе, канд. физ.-мат. наук Алесе Сергеевне Ибрагимовой за консультации при реализации алгоритмов и ценные советы в ходе исследований, а также Станиславу Андреевичу Боярскому и Роману Евгеньевичу Семенко за обсуждения и подцержку при выполнении работы.

СОДЕРЖАНИЕ РАБОТЫ

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

2

+ СКУ(Л) =0, Л4 + -Ча = Я<Э + сп Л + Сх21,

О

<7( + <Ну(1) = (Л, <3) + сР, 14 + У(аг К) = + С21Л + с221.

о

Здесь К - электронная плотность, Е - энергия электрона, Л = Ди, I = Яq, и = (и(х), и^)- вектор скорости электронов (за основу принята декартова

система координат (х,у)), q = - поток энергии, Р = —

а — НЕ, аз = <3 = = ((рх,<ру), <р = х, у) - электрический

потенциал, удовлетворяющий уравнению Пуассона:

= <рхх + <руу =/3(Л-р); (2)

р = р(ж, у) - плотность легирования (заданная во внутренней области полупроводника функция). Коэффициенты сц, ..., С22, с системы (1) являются гладкими функциями от энергии Е, Р > 0 - постоянная. В выкладках мы используем также диэлектрическую постоянную е = Существенной сложностью при поиске численных решений (1), (2) является наличие в модели малых параметров и функций, таких как е и р(х,у).

Значительная часть введения посвящена формулировке основных целей работы и идей, служащих фундаментом разработанных методов, обоснованию актуальности указанных целей и идей, изложению содержания работы и обзору полученных результатов - приведено схематическое описание двух новых подходов для поиска стационарных решений гидродинамической МЕР модели и созданной технологии построения вычислительных моделей для поиска решений Ш и 2Б нестационарных смешанных краевых задач.

Первая глава диссертации состоит из шести параграфов. Здесь на примере модельной задачи описан и обоснован один способ поиска стационарных решений уравнений гидродинамической МЕР модели (1), (2). Этот способ (в работе мы именуем его «первый подход», см. рис. 1) основан на сведении соотношений (1), (2) в стационарном случае к системе из трёх уравнений Пуассона и применении к ним нестационарной регуляризации. При построении вычислительной модели данного подхода мы разработали

(!)

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

В §1 система уравнений МЕР модели (1), (2) в стационарном случае сведена к трём уравнениям Пуассона. В §2 для уравнения Пуассона

Л*,УФ = Ф** + ФУУ = f{x,y), (а:,у) € О (3)

рассматривается модельная краевая задача Неймана-Дирихле в прямоугольной области. Здесь Ф(х, у) - неизвестная функция, f(x, у) - известная правая часть. Также в §2 получены априорные оценки на норму Ф(х,у).

Для поиска решений модельной задачи в п. 3.1-3.3 §3 предлагаются три вида нестационарной регуляризации:

параболическая регуляризация щ = Ли — f(x, у): регуляризация Соболева щ — Лщ = Ли — f(x,y); гиперболическая регуляризация utt + Кщ — Ли + f(x, у) = 0. Здесь u(t, х, у) - новая неизвестная функция, К > 1 - постоянная, t -временная переменная. Целесообразность подобного перехода от уравнения Пуассона к нестационарному регуляризованному уравнению оправдана, потому что это позволило нам реализовать идею о поиске стационарного решения уравнений МЕР модели с помощью метода установления. При этом стационарное решение ищется как предел u(t, х, у) при t —^ оо.

Обоснование сходимости метода установления в случае применения каждой из трёх регуляризаций представляет отдельную задачу, так как прежде всего необходимо доказать однозначную разрешимость регуляризованного уравнения, а затем асимптотическую устойчивость по Ляпунову его решения при t —► оо. С этой целью в п. 3.1-3.3 §3 для каждого из трёх видов регуляризаций получены априорные оценки на нормы решений регуляризо-ванных задач, из которых следует сходимость u(t, х, у) —> Ф(х, у) при t —> оо, а также однозначная разрешимость и устойчивость стационарного решения.

В §4 гл. 1 для поиска приближённых решений регуляризованных уравнений предложена дифференциально-разностная модель, основанная на методе прямых и принципах схем без насыщения. При этом ихх аппроксимирована второй производной от интерполяционного полинома

Пх,и) = I Eyr\JllSXi^Nx)ui{t,yb 0 <х <п

с узлами Xi, i = 1, ...,N в нулях многочлена Чебышева.

В §4 проведена дискретизация регуляризованных уравнений по переменной t. Вводя обозначения: и" (у) = и(пЛ, у) = и (у), п = 0,1,..., и(у) = = ип+1(у), Д - шаг разностной сетки по времени, и' = щ, в итоге приходим к системе ОДУ второго порядка

Y" = BY + F, (4)

Y = ({¿i(t, y),..., uw(t, у))т. Выражения для элементов матрицы В и компонент вектора Т зависят от вида применённой регуляризации.

К системе (4) добавляем граничные условия, которые следуют из условий Неймана-Дирихле модельной задачи §2.

Важно отметить, что при использовании параболической регуляризации для доказательства устойчивости дифференциально-разностной модели из ОДУ в §4 гл. 1 нам удалось вывести разностный аналог априорной оценки, полученной для функции Ф в §2. Таким образом, по крайней мере, в случае модельной линейной задачи для уравнения (3) мы обосновали адекватность созданной вычислительной модели по отношению к рассмотренной задаче.

Параграф 5 гл. 1 посвящён поиску приближённых решений системы (4), с граничными условиями в виде интерполяционного кубического сплайна класса С2:

К2

S(у) = (1 - r)Yfe + rYfc+i - -f?(l - т)[(2 - т)тк + (1 + rjmfc+r], (5)

о

где

~ _ У Ук ^ у € [yfci ^+1]( yk = i(fl к = 0,К — 1, Khy = пу о

Yk=Y(yk), mk=Y(yk).

С учётом (5) из (4) была получена трёхточечная матричная схема

{/дг - - + ^Б}>Yfc + |Гдг - f fi}Yfc+1

h2 _

= -^{^•-1+4^+^+1}, k = l,K-l,

Y0=Y1,YK = LYK_1, (7)

где = Ik), L ~ диагональная матрица граничных условий порядка N. Решение (6), (7) на каждом временном слое несложно найти с помощью хорошо известного метода матричной прогонки.

Алгоритм, построенный в гл. 1, был запрограммирован на языке Delphi 6 (среда Object Pascal). В §6 продемонстрировала эффективность предложенной вычислительной схемы при решении тестовой краевой задачи для уравнения Пуассона с известным точным решением. При этом исследована зависимость относительной погрешности численного решения тестовой зги-дачи от различных параметров.

Вторая глава диссертации содержит три параграфа. Она посвящена построению и обоснованию вычислительной модели подхода 2 (см. рис. 1).

(6)

С этой целью система (1) преобразована к симметрической форме и для поиска её решений сконструирован класс «устойчивых» разностных схем.

В §1 гл. 2 система уравнений (1) на её гладких решениях записана в недивергентной форме:

Здесь U =

Ut + ßUx + CUy = F.

{ 0 \ RQ + ciiJ + c12I

(J.Q) + cP

\ fffQ + C21J + C22I / В, С - матрицы определённого вида.

(8)

F =

( JW \

• J = { j(») J = fíu'

Определение. Вещественная матрица А = А* = (ау), г, 7 = 1,6 называется симметризатором для системы (8), если:

1 )А>0, 2)В = АВ = В\ 3) С = АС = С*.

Здесь Л* - транспонированная матрица А.

В §1 для системы (8) построен симметризатор, после чего соотношение (8) преобразовано к симметрической ¿-гиперболической по Фридрихсу форме:

А • и£ + В • их + С • иу = АР. (9)

Второй параграф гл. 2 посвящён выводу априорной оценки на норму решения системы (8) в предположении, что электрический потенциал <р{Ь, х) - известная достаточно гладкая ограниченная функция. Прежде всего в §2 при некоторых допущениях был получен другой вид соотношения (9):

(AU)t + (BU)X + (CU), = +

(10)

где V — V + С, V = АО, О - матрица коэффициентов сц,...,г.22, с, С — Вх + Ву, Т - вектор определённого вида.

Далее, умножая системы (9), (10) скалярно на вектор и и складывая полученные выражения, мы в итоге нашли:

(U, AU)t + (U,BU)X + (U,CU)„ = (U, AU) + 2(и,Я-

(И)

Здесь Л = V + V* 4- (С + С)/2.

Считая, что система (1) имеет ограниченное гладкое решение U(f, х, у) в области П = {(£, x,?/)|To<i<T< 00, (х, у) £ В2}, где Т0 > 0 - достаточно большое число, полагая:

1) |U|2 = (U,U) 0 при у/х2 + у2 -» оо;

2) нормы матриц В, С ограничены в области П;

3) в области П матрица А > 0, R > 0, Е > 0;

4) в области П вектор-функции Q, u, q и функция Е - ограниченные и гладкие,

и интегрируя (11) по В.2, мы вывели априорную оценку на норму решения системы (1):

Щ < 1{Т0)ехр{Д(4 - Го)}, Г0 < 4 < Г < оо. (12)

Здесь Д - положительное число, зависящее от спектра матриц Л и Л и от величины ае[(3[, 1(Ь) = / (и, А\5)йхйу - интегральная норма решения. в?

В §3 гл. 2 для поиска решений системы (8) в области П предложен один класс «устойчивых» разностных схем. Идея, заложенная нами в основу созданного класса, отражает принцип адекватности и состоит в том, чтобы получить для разностных схем аналоги выражений (11), (12), доказанных дчя дифференциальной задачи в §2.

Для записи разностных отношений в области П построена сетка с шагами Д по времени, кх, ку по переменным х, у. Введены такие обозначения:

= и(пД, аЬ) = иа = ив, = иау = и - сеточная вектор-функция;

а = (ах,ау) - мультииндекс, Ь = (/г1,/гу),аЬ = (ахЬх,ауЬу), Ы, Ы = 0,1, ...;п = N0,..., ЛГ; А = ДЛ/о,Г = ДЯ1;

X, Ф®, Фу, Ф*1, Фу1 - операторы сдвига: = и^"1 = О, Ф^и =

= Ф?1 и = = Фя,у);

г, - разностные операторы: т — х~ 1> £с = Ф* - 1>

& = Фв - 1.1« = 1 - Ф*\ !у = 1- 9-1.

Как известно, симметрические матрицы В, С можно представить в виде В — В+ - В-, С = С+ — С_, где В±,С± > 0 - некоторые симметрические матрицы. В §3 предложен конструктивный способ поиска матриц В±,С±:

В± = й(Е)А ± В/2, С± = ¿(Я) А ± С/2.

Здесь ¿(Е) > / —- некоторая функция, выбираемая так, чтобы у

обеспечить выполнение неравенств В±,С± > 0.

Предлагаемый класс конечно-разностных схем имеет следующий вид:

+ г, {К'В+(Е^хи + К^ив^щ - КВ-{£Ы)Ьи--ка1+1ав-(£1хЩ} +гу {кс+(Е^%и + кау^у(с+(Е^Щ--КС-(£Ь%О - £в,+1£р(С_(£&))0)} = АК(2Э + £)© + 2АК&,

где

А = А(Е), А = A(Ê), V = V(Ê), £ = £(Û), Т =

К, К, Kat±u Кау± 1 - диагональные матрицы неизвестных;

тх — у = Е(х\,Eiyï,- «промежуточные» значения

энергии fco > 0- некоторое число,

Замечание 1. Конечно-разностное выражение (13) определяет целый класс конечно-разностных схем, поскольку возможны самые различные способы определения «промежуточных» значений энергии Е.

В §3 для сеточных функций и матриц сформулированы аналоги предположений 1) — 4) §2. При использований этих предположений после домно-жения (13) на вектор V = (1,1,1,1,1,1)т скалярно получен разностный аналог равенства (11), из которого при подходящем значении ко, в условии малости шага Д выведен аналог оценки (12):

/ 1 \ п-n0

In0,ti = No + 1,...,N. (14)

Здесь £ - положительное число, зависящее от спектра матриц А(Е£), A(U£) и величины as|Q2|, In - норма сеточной вектор-функции UJj определённого вида.

Замечание 2. Из (14) следует также «устойчивость» разностной модели (13) в энергетической норме \/7„. Мы взяли слово «устойчивость» в кавычки, ибо, строго говоря, ещё требуется доказать, что на решениях конечно-разностной вычислительной модели (13) выполнены условия 1) — 4) §2 и их разностные аналоги.

Третья глава диссертации включает шесть параграфов. Здесь для уравнений МЕР модели в 1D случае сформулирована хорошо известная из физики полупроводников задача о баллистическом диоде. При поиске стационарных решений данной задачи продемонстрирована высокая эффективность подходов 1 и 2 (см. рис. 1), детально описанных в предыдущих главах. В этой главе мы также предложили несколько других способов поиска приближённых решений задачи о баллистическом диоде.

В §1 рассмотрены уравнения МЕР модели в одномерном случае в безразмерном виде. Для них дана постановка задачи о баллистическом диоде (рис. 2), т.е. заданы условия на границах г = 0иг = 1и начальные данные при t = 0. Характерными величинами задачи являются постоянные L и N+, где L - длина диода, N+ - плотность легирования п+-зоны (см. рис. 2). При этом безразмерная плотность легирования р(х) такова:

!1, если х принадлежит п+-зоне,

г N

о = тгг, если х принадлежит п-зоне,

Рис. 2. Схематическое представление п+ — п — п+ баллистического диода

где /V - плотность легирования п-области. Типичный сглаженный профиль функции р(х) изображен на рис. 2.

Параграф 2 состоит из четырёх пунктов. В нём при поиске стационарных решений задачи о баллистическом диоде мы следовали идеям гл. 1 (см. рис. 1, подход 1). В п. 2.1 по аналогии с 2В случаем стационарные уравнения Ш МЕР модели сведены к системе из трёх уравнений второго порядка. В п. 2.2 рассмотрена гиперболическая регуляризация полученных уравнений. Важно отметить существенные отличия исследований, проведённых нами для Ш и 2Б задач. Если в 2Б случае все обоснования алгоритма были проведены на примере модельной линейной задачи для уравнения (3), то в п. 2.2 §2 (при некоторых ограничениях на функцию р(х) и постоянную Ъ) мы вывели априорную оценку на норму решения регуляризованной задачи на исходном квазилинейном уровне. Как и ранее, эта оценка обосновывает применимость метода установления и доказывает асимптотическую устойчивость стационарного решения. В п. 2.3 мы рассмотрели модельную задачу для регуляризованных уравнений. Основываясь на технологии, описанной в гл. 1, мы свели её к поиску решения ОДУ второго порядка с краевыми условиями, которое записали в виде интерполяционного кубического сплайна класса С2. В результате был получен вычислительный алгоритм для поиска решений регуляризованных уравнений с помощью метода матричной прогонки. В п. 2.4. на основе данного алгоритма и метода установления построена схема поиска стационарных решений задачи о баллистическом диоде.

Параграф 3 гл. 3, включающий три пункта, посвящён поиску решений регуляризованных уравнений способами, отличными от предложенных в §2. С этой целью в §3 разработаны вычислительные алгоритмы, основанные на методе ортогональной прогонки (п. 3.1), технике интегральных соотношений (п. 3.2) и схеме предиктор-корректор (п. 3.3).

Параграф 4 гл. 3 состоит из двух пунктов. В данном параграфе подход 2 (см. рис. 1) реализован в Ш случае. Важно отметить, что ход рассуждений и выкладки, необходимые для создания класса «устойчивых» разностных схем подхода 2, перенесены с 2В на Ш случай с минимальными изменениями (основная идея же осталась неизменной, см. п. 4.1 §4).

Далее, в п. 4.2 описан вычислительный алгоритм для поиска стацио-

нарных решений задачи о баллистическом диоде с помощью метода установления. При этом класс «устойчивых» разностных схем для 1D задачи преобразован к трёхточечной матричной схеме. В соответствии с задачей о баллистическом диоде для неё поставлены граничные условия и решение на каждом временном слое найдено методом матричной прогонки.

В §5 гл. 3 с целью верификации методов, предложенных в §2-4 гл. 3, рассмотрена задача о баллистическом диоде при £ = 0 (е - диэлектрическая постоянная).

В §6 гл. 3, состоящем из пяти пунктов, речь идёт о создании на базе алгоритмов, описанных в главе, прикладных программ. Важно отметить, что все разработанные в гл. 3 алгоритмы удалось реализовать на ЭВМ. При этом на языках Delphi 6 (среда Object Pascal) и Java были написаны программы для нахождения приближённых решений задачи о баллистическом диоде. В п. 6.1 приведены значения параметров и переменных задачи, необходимые для расчётов. В п. 6.2 подробно описаны проблемы, возникшие при вычислениях, обусловленные наличием малых параметров, больших градиентов и нелинейностью уравнений модели. Для решения этих сложностей в п. 6.2 предложено использовать итерации tío нелинейности и нелинейное сглаживание. В п. 6.3 мы остановились на алгоритме, использующем класс «устойчивых» разностных схем, и описали эффективные средства, благодаря которым удалось получить решение с необходимой точностью. Эти средства включают применение прогонок в разных направлениях и расчёт потенциала двумя способами - с помощью метода прогонки и с помощью квадратурной формулы трапеций. Алгоритм, основанный на «устойчивых» разностных схемах, удалось распараллелить, после чего было проведено несколько быстрых расчётов на кластере НГУ (см. https://www.nusc.ru/). Детали распараллеливания, а также результаты вычислений расположены в п. 6.4. В п. 6.5 представлены графики численных решений, полученные при расчётах. Здесь же проведено сравнение результатов и эффективности предложенных в гл. 3 алгоритмов. На рис. 3 изображены графики энергии электронов Е, электрического потенциала tp, плотности электронов R и скорости электронов и, полученные в ходе расчётов с использованием алгоритма подхода 1 и схемы подхода 2 (см. рис. 3 а, б).

Четвёртая глава состоит из четырёх параграфов. Здесь речь идёт о применении технологии, разработанной в гл. 1, для поиска стационарных решений задачи о переносе заряда в 2D кремниевом транзисторе MOSFET с наноканалом из оксида кремния. Подробное описание такого транзистора с электронной проводимостью приведено в работе V. Romano6.

В §1 для уравнений (1), (2) дана постановка задачи о переносе заряда в транзисторе MOSFET. Схематическое изображение данного транзистора

6Romano V. 2D Numerical Simulation of the МЕР Energy-Transport Model with a Finite Difference Scheme // J. Comp. Fhys. 2007. Vol. 221. P. 439-468.

Рис. 3. Результаты расчётов, полученные: а - с применением подхода 1 при V = 1 Volt, 5 = 0,004, L — 6 х 10~7m; б - при использовании схемы подхода 2 для V = 1,5 Volt, 5 = 0,004, L = 3 х 10~7т

приведено на рис. 4. Его существенной особенностью является наличие на-ноканала из оксида кремния О. с-

Замечание 3. Поскольку перенос заряда в наноканале Qg (см. рис. 4) отсутствует, то в области Яс электрический потенциал Ф(4, х, у) удовлетворяет уравнению Далласа:

у

ííb(sio2) (т;

р.i ÍÍ+

£ р

p-ô fi (Sil

1/4 5,1« ЬШк Н/16 3/4

ДФ = Фхх + Фуу = 0, (15)

В §1 гл. 4 для математической модели (1), (2), (15) выставлены граничные условия. Следуя схеме рассуждений гл. 1, мы свели поставленную краевую задачу к поиску решения системы трёх уравнений Пуассона для неизвестных Д, <р, $ =\Е — 1в области Г2 = {(х, у) : 0 < х,у < 1}.

В §2 гл. 4 предложен один способ преодоления сложностей, возникших на пути применения технологии гл. 1, для поиска приближенных решений полученной краевой задачи для

Рис. 4. Схематическое представление 2D кремниевого транзистора MOSFET

уравнений Пуассона. Непосредственному применению данной технологии мешало только одно обстоятельство - условия склейки на границе S -|ФУ = tpy, Ф = (р, связывающие потенциал <р с потенциалом Ф. Пусть 1Х и 1У

- длина и ширина наноканала ÎIg соответственно. В §2 показано, что если отношение Ем = 'f- пренебрежимо мало, то краевую задачу для потенциала ср в области SÏ можно доопределить, т. е. на множестве S можно сформулировать дополнительное смешанное краевое условие для функции ¡р. Вид этого условия (см. (16)) определяется в результате упрощения процедуры поиска потенциала Ф в наноканале Па-

ф, 1) + 31у<ру(х, 1) = G, (х, 1) € S. (16)

В результате, учитывая (16), мы можем использовать технологию гл. 1 для поиска приближённых значений потенциала ip в области il

В §3 гл. 4 для верификации технологии, предложенной в гл. 1, в совокупности с методом установления и параболической регуляризацией мы применили хорошо известный метод продольно-поперечной прогонки (п.п.п.). Этот метод основан на использовании разностных отношений, посредством которых аппроксимируются производные неизвестных функций. Для его реализации необходимо ввести пространственную и временную разностные сетки. Переход от предыдущего временного слоя к следующему осуществляется в два этапа: первый этап - продольная прогонка, второй

- поперечная. В соответствии с идеей метода установления, данные операции осуществляются до тех пор, пока решение не установится. В результате находим стационарное решение задачи о переносе заряда в MOSFET.

Параграф 4 посвящён проблемам реализации построенных вычислительных схем. Здесь описано несколько высокоэффективных способов, которые позволили нам добиться сходимости метода установления и получить решение задачи о переносе заряда в 2D транзисторе MOSFET с необходимой точностью. Эти способы включают итерации по нелинейности, фильтр нелинейного сглаживания, применённый вдоль осей абсцисс и ординат, и процедуру «вытягивания параметров».

В §4 приведены результаты расчётов, проведённых с использованием разработанной технологии и метода п.п.п. После трансформации вычислительной схемы подхода 1 с использованием новых вспомогательных переменных были найдены стационарные решения задачи о переносе заряда в транзисторе MOSFET (рис. 5). Данные решения были получены при следующих значениях параметров:

VD = IV, VG = IV, В = -25,328, <5 = -0,001, ly = 20пт.

Здесь Vd - напряжение на стоке drain, Vg - напряжение на затворе gate, В - безразмерное напряжение на корпусе bulk, S - безразмерная плотность легирования в области Q \ Î2+ (см. рис. 4).

Необходимо заметить, что получить решение для данного набора параметров задачи удалось при помощи метода п.п.п., хотя при расчётах с дру-

Рис. 5. Распределения энергии и электрического потенциала в транзисторе МОЭРЕТ, полученные в ходе расчётов методом п.п.п.

гими значениями параметров алгоритм, основанный на идеях гл. 1, оказался намного эффективнее прогонки.

В Приложении I диссертации подробно описан процесс обезразмерива-ния уравнений МЕР модели, приводящий к системе (1), (2). Приложение II содержит явные выражения для коэффициентов с, сц, ...,С22, стоящих в правой части системы (1).

В Приложении III предложенная в диссертации технология апробирована при изучении вопроса об устойчивости ударных волн в сжимаемом вязком газе, движение которого описывается известными уравнениями Навье-Стокса. Приложение III содержит 3 параграфа, в которых дана постановка задачи для уравнений Навье-Стокса и сделаны необходимые преобразования. С помощью разработанной технологии были организованы расчёты, результаты которых (см. §2, 3 Приложения III) показали устойчивость рассмотренных ударных волн при определённых значениях параметров.

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

В §1 Приложения IV кратко описана система гидродинамических уравнений для газосодержащих нефтяных слоистых систем. С учётом весьма правдоподобных упрощающих предположений после преобразований данной системы для малых деформаций и слоистой структуры поставлена двумерная нестационарная смешанная задача.

Для поиска численного решения полученной задачи в §2 использована технология, разработанная в диссертации. В результате при определённых значениях параметров обнаружен эффект роста со временем амплитуды ко-

лебания решения (рис. 6). Это и есть искомый параметрический резонанс, приводящий к разрушению слоистой структуры.

10 000 8 000 6 ООО 4 000 2 000

-2 000 •4 000 -6 000 -в 000 •10 000

О 1847,100000 4282,170000 6717,240000 9152,31000 11280

Рис. 6. Зависимость амплитуды колебаний решения от времени

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

ПУБЛИКАЦИИ ПО ТЕМЕ ДИССЕРТАЦИИ

В рецензируемых журналах, рекомендуемых ВАК:

1. Блохин А. М., Ибрагимова А. С., Семисалов В. В. Конструирование вычислительного алгоритма для системы моментных уравнений, описывающих перенос заряда в полупроводниках // Мат. моделирование. 2009. Т. 21. N 4. С. 15-34.

2. Блохин А. М., Ибрагимова А. С., Семисалов Б. В. Конструирование вычислительных алгоритмов для задачи о баллистическом диоде // Журн. выя. мат. и мат. физ. 2010. Т.50, N 1. С.188-208.

3. Блохин А. М., Боярский С. А., Семисалов Б. В. Один способ построения разностных моделей для системы моментных уравнений, описывающих перенос зарядов в полупроводниках // Вестн. НГУ. Серия: Математика, механика, информатика. 2009. Т.9, вып. 4. С. 3-15.

4. Блохин А. М., Семисалов Б. В. Конструирование одного класса вычислительных схем в задаче о баллистическом диоде // Мат. моделирование. 2010. Т. 22(5). С. 3-21.

В зарубежных журналах:

5. Blokhin А. М., Boyarsky S. A., Semisalov В. V. On an approach to the construction of difference schemes for the moment equations of charge transport in semiconductors // Le matematiche. 2009. Vol. LXIV. Fase. I. P. 77-91.

6. Blokhin A. M., Semisalov В. V. Design of numerical algorithms for the problem of charge transport in a 2D silicon MOSFBT transistor with a silicon oxide nanochannel // J. Phys.: Conf. Ser. 2011. Vol. 291. Art. 012016.

URL: http://iopscience.iop.Org/1742-6596/291/l/012016

В трудах международных конференций и семинаров:

7. Blokhin А. М., Semisalov В. V. The construction of one class of numerical algorithms in ballistic diode problem // Int. School and Seminar on Modern Problems of Nanoelectronics, Micro- and Nanosystem Technologies Proc. «INTERNANO-2009». NSTU, Novosibirsk, Russia - Oct. 28-31.

8. Blokhin A.M., Semisalov В. V. The construction of numerical algorithms for the ballistic diode problem // Proc. Int. Conf. on Computational Technol. In Electrical and Electronics Engineering «Sibircon-2010». Irkutsk, Listvyanka, July 11-15, 2010. Vol. I. P. 432-437 (IEEE Region 8).

9. Блохин A. M., Семисалов Б. В. Конструирование вычислительных алгоритмов в задаче о переносе заряда в кремниевом транзисторе MOSFET с нанока-налом из оксида кремния: Тез. докл. III Междунар. форума по нанотехнологиям. М.,1-3 нояб. 2010 (электронный ресурс).

В тезисах конференций:

10. Ибрагимова А. С., Семисалов Б. В. Конструирование вычислительного алгоритма для одномерной МЕР-модели переноса зарядов в полупроводниках // XLVII Международная студенческая конференция (ISSC). Россия. Новосибирск, 12.04 - 15.04.2009. С. 222.

11. Семисалов В. В., Боярский С. А. Один способ построения разностных моделей для системы моментных уравнений, описывающий перенос заряда в полупроводниках // XLVIII Международная студенческая конференция (ISSC). Россия. Новосибирск, 10.04 - 14.04.2010. С. 220.

Научное издание

Семисалов Борис Владимирович

Математическое моделирование в задачах переноса заряда в полупроводниковых кремниевых устройствах

05.13.18 - математическое моделирование численные методы и комплексы программ

Автореферат диссертации на соискание учёной степени кандидата физико-математический наук

Подписано в печать 05.05.2011 г. Формат бумаги 60x84 1/16. Уч.-изд. л. 1,3. Усл. печ. л. 1,3. Тираж 100 экз. Заказ N0.139 Редакционно-издательский центр НГУ. 630090, Новосибирск-90, ул. Пирогова, 2.

Введение

Глава 1. Регуляризованная дифференциально-разностная модель (Подход 1)

§1. Преобразование уравнений МЕР модели к системе уравнений

Пуассона.

§2. Модельная краевая задача для уравнения Пуассона.

§3. Нестационарные регуляризации.

3.1. Параболическая регуляризация.

3.2. Регуляризация Соболева.

3.3. Гиперболическая регуляризация.

§4. Устойчивая дифференциально-разностная модель.

§5. Вычислительный алгоритм для модельной задачи.

§6. Реализация алгоритма для тестовой задачи.

Глава 2. Конструирование класса «устойчивых» разностных схем (Подход 2)

§1. Сведение уравнений МЕР модели к симметрической по

Фридрихсу системе.

§2. Априорная оценка для системы (1.6).

§3. Один класс «устойчивых» разностных схем.

Глава 3. Задача о баллистическом диоде

§1. Постановка задачи.

§2. Применение подхода 1 в Ш случае.

2.1. Об аналогах с 2В случаем.

2.2. Глобальная априорная оценка для регуляризованной задачи.

2.3. Вычислительная модель подхода 1 в Ш случае.

2.4. Схема поиска решения задачи (2.12)-(2.15), (2.11).

§3. Другие методы поиска решения модельной краевой задачи (2.28), (2.29).

3.1. Метод ортогональной прогонки.

3.2. Преобразование к системе интегральных уравнений.

3.3. Схема предиктор-корректор.

§4. Применение подхода 2 в Ш случае.

4.1. Об аналогах с 2Т) случаем.

4.2. Описание вычислительного алгоритма.

§5. Стационарная задача в случае г = 0.

§6. Реализация вычислительных алгоритмов.

6.1. Значения параметров и переменных, необходимые для расчётов.

6.2. Итерации по нелинейности. Нелинейное сглаживание.

6.3. Тонкости реализации алгоритма (4.12'), (4.14).

6.4. Распараллеливание алгоритма (4.12'), (4.14).

6.5. Обсуждение и.сравнение результатов.

Глава 4. Задача о переносе заряда в 2D транзисторе

MOSFET

§1. Постановка задачи.

§2. Дополнительное краевое условие для электрического потенциала.

§3. Метод продольно-поперечной прогонки.

§4. Реализация вычислительных алгоритмов.

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

Моделирование процесса переноса зарядов в полупроводниковых приборах субмикронного размера основывается на кинетическом уравнении Больцмана, описывающем движение носителей зарядов (электронов и дырок). Для функции распределения электронов / = /(£, х, V) в трёхмерном случае это уравнение имеет вид

Здесь х = (х1,х2,хз), V = (г^г^^з) - векторы координат и скорости, д - заряд электрона, т* - эффективная масса электрона, Е = {Е\, Е2, -вектор напряжённое™ электрического поля, С} - оператор Больцмана, учитывающий взаимодействие электронов с атомной решёткой. Уравнение (1) написано для случая, когда имеется только электронная проводимость. При наличии нескольких типов носителей зарядов (электронов и дырок, «легких» и «тяжелых» дырок и т.д.) необходимо ввести свою функцию распределения для каждого типа. В результате мы будем иметь столько кинетических уравнений вида (1), сколько есть таких носителей. При этом в правой части формулы (1) появится сумма по всем типам частиц, на которых рассматриваемые носители заряда могут рассеиваться. Таким образом, если имеет место взаимное рассеяние или превращение частиц разных типов, то мы приходим к системе связанных кинетических уравнений.

Для поиска приближённых решений уравнения (1) предлагается множество различных способов. Наиболее распространённый из них - метод Монте-Карло. Однако этот метод обладает рядом недостатков: он требует существенных вычислительных затрат и в некоторых случаях (например, если концентрация носителей заряда в отдельных областях полупроводникового устройства мала) не позволяет получить точный результат. Поэтому возникает потребность в других моделях, представляющих собой разумный компромисс между физической точностью и вычислительной эффективностью. Например, приемлемая точность может быть достигнута при решении уравнений переноса, полученных для моментов уравнения Больцмана, таких как п(£,х) = У м3

7Ш(г,х) = J г м2 гае(*,х) = J и|2

V — 71

Здесь п - концентрация, и - средняя скорость, е - внутренняя энергия носителей зарядов, \\\2 = (V, V) и т.д.

Простейшей моделью переноса заряда, которая получена методом моментов из уравнения Больцмана, является дрейф-диффузионная модель. Классическая дрейф-диффузионная модель физики полупроводников была предложена в 1949-1950 гг. ЭсЬосЫеу [1] и уап КооэЬгоеск'ом [2]. Она состоит из уравнения Пуассона для электрического потенциала: --(р - п + ^ - дга) (2) и уравнений неразрывности для носителей заряда: дть др - сНу.^ = Я(р, гс), — - сМр = Я(р, п), (3) где п(£, х), х) - концентрации электронов и дырок соответственно. Векторы плотности электронного и дырочного потенциалов записываются в таком виде:

Зп = ИгУп - 1лппЧ(р, Зр = ИрУр - ЦррЧф.

Здесь Ип, ]Эр - коэффициенты диффузии, /1Р - подвижности электронов и дырок соответственно, во ~ относительная диэлектрическая проводимость полупроводника. Предполагается, что материал полупроводника легирован донорной и акцепторной примесями с концентрациями Л^(х), ]Уа(х) и в нём происходит рекомбинация частиц со скоростью Н(р,п).

Необходимо отметить, что дрейф-диффузионная модель достаточно полно изучена в математическом смысле и сегодня используется в основной массе программ, предназначенных для моделирования процесса переноса заряда в полупроводниковых устройствах. Математические аспекты модели (2), (3) исследованы многими авторами. Отметим в этой связи работы [3]-[5].

М.Э. Моек, полагая в (2), (3) £>п = = рп = рр = 1 при некоторых ограничениях на область определения и начальные данные, доказал существование единственного решения задачи Неймана для уравнений (2), (3), гладким образом зависящего от начальных данных (см. [6]). В работе [7] проведён анализ стационарной смешанной задачи Неймана-Дирихле для (2), (3). Вопросы точности описания дрейф-диффузионными моделями различных физических явлений обсуждаются в [8]. В статье [9] получены обобщённые уравнения диффузии и дрейфа носителей заряда в неупорядоченных полупроводниках. Главная особенность этих уравнений - наличие частной производной по времени дробного порядка, что позволяет описывать нормальный и дисперсионный перенос в рамках единой математической модели. Численное же моделирование процесса переноса заряда в полупроводниковых устройствах восходит к известной работе [10], авторы которой предложили устойчивую дискретизацию дрейф-диффузионных уравнений, используемую и по сей день.

Однако возрастающая миниатюризация современных электронных приборов требует более точного описания процесса переноса энергии в полупроводниках, что имеет первостепенную важность для описания таких феноменов, как «горячие» электроны, ударная ионизация, генерация тепла и т.п. Следовательно, возникает необходимость расширения общепринятой дрейф-диффузионной модели, с целью учёта энергии носителей заряда. Эта цель достигается в так называемых гидродинамических моделях.

При построении гидродинамических моделей выбирается подходящая процедура замыкания, которая позволяет из бесконечной системы моментов уравнения Больцмана получить замкнутую систему из конечного числа уравнений. Многообразие процедур замыкания и определяет весь спектр гидродинамических моделей (обзор некоторых из них приведён в [11]-[13]).

Одна из самых первых гидродинамических моделей получена Вкйек^аег'ом (см. работу [14]) и изучалась Вассагаш, Уогс1ета11!ом (см. [15]) и другими авторами. При наличии одного типа носителей зарядов (электронов) модель записывается в виде системы для электронной плотности р, скорости и и плотности энергии 8 (см. [12]): + div(pu) = 0, (4) + div(pu <g> u) + Vp(p, Т) = pVV - ри, (5) д£ + div(£u + р(р, Т)и - ^VT) = ри • W - (S - (6)

AV = p-C(a?). (7)

Здесь Т - температура, р(р,Т) - давление, V - электрический потенциал, С(х) - профиль легирования, к в ~ постоянная Больцмана, ТУ, - температура решетки, = |рквТь• Плотность энергии имеет следующее выражение: где m - эффективная электронная масса.

Однако полный математический анализ гидродинамической модели (4)-(7) пока отсутствует. В то же время хорошо изучена упрощенная, так называемая изотермическая гидродинамическая модель, которая получается из (4)-(7) в предположении, что температура постоянна: + div(pu) - 0, (8) div(pu <8> u) + Vp = pW - ри, (9)

С/6

AV = p-C(x), (10) где р — р(р). Часто предполагается, что р(р) = ^р7, 7 > 1.

Для изотермической модели в работах [16]-[18] исследовались такие вопросы, как существование и асимптотическое поведение решения. Также имеется большое количество работ, посвященных вычислительным аспектам модели (см., например, [19]). Некоторые численные методы для поиска приближённых решений (4)-(7) в стационарном случае предложены в статье [20]. Работа [21] посвящена моделированию ударных волн в стационарном случае для субмикронных полупроводниковых устройств. В статье [22] осуществлён сравнительный анализ изотермической модели (8)—(10) и дрейф-диффузионной модели в одномерном случае, в рамках которого изучена асимптотика решений сравниваемых моделей при t —>■ оо. Также в [22] для упомянутых моделей исследованы вопросы существования и единственности решения задачи Коши и начально-краевой задачи.

В настоящей диссертации рассматривается одна из последних моделей гидродинамического типа, которая была предложена итальянскими физиками Anile и Romano [23, 24]. Авторы модели отказываются от предположения идеальной среды (когда тензор напряжения считают шаровым) и от обычно применяемого для замыкания гидродинамических моделей закона Фурье для потока тепла:

Q = -kVT

Модель Anile и Romano представляет собой квазилинейную систему уравнений, записанных в форме законов сохранения. Эти законы получены из системы моментных соотношений для уравнения переноса Больцмана с помощью так называемого принципа максимума энтропии (или МЕР от Maximum Entropy Principle). Модель имеет следующий вид: + div(nV) = 0, д{пР) dt dnW dt d(nS) dt V -nW + neE = nCp, div(nS) + ne( E, V) = nCw,

10 WT 5 öne~ГЕ = nCW

6 m*

П) и рассматривается совместно с уравнением Пуассона бДФ = e(n - N).

12)

Здесь п, V, И7, Э - соответственно электронная плотность, средняя скорость электрона, средняя энергия электрона, поток энергии; Р = т*У - средний момент кристалла, е - абсолютное значение заряда электрона, Е = —УФ -электрическое поле, Ср(И^), Сцг{УУ), С\у(^0 - члены производства балансных уравнений, е - диэлектрическая постоянная, А^», Ад - плотности доноров и акцепторов, N = N0 — А/д.

Замечание 1. В данной работе нам потребуются уравнения (11), (12) в безразмерном виде (процесс обезразмеривания (11), (12) подробно описан в Приложении I и в [25]):

Щ + сПУ(Л) = О, 2 - У(7 = + СцЗ + С121, О а!у(1) = (л, <э) + ср, 5 У(ееЛ) = -сгС^ + с21Л + с221. О

13)

Здесь Н - электронная плотность, Е - энергия электрона, Л = Яи, I = Яц, и = - вектор скорости электронов (за основу принята декартова система координат (х,у)), д = ~ поток энергии, Р = Е — ,

7 = ЯЕ, аз = ЩЕ2, = У^? = (срх,(ру), <р - - электрический потенциал, удовлетворяющий уравнению Пуассона: фхх + <Руу = - р)\ (14) р = р(х, у) - плотность легирования (заданная во внутренней области полупроводника функция). Коэффициенты сц, ., С22, с системы (13) являются гладкими функциями от энергии Е, выражения для которых приведены в Приложении II, ¡3 > 0 - постоянная (по поводу постоянной (3 см. [25]). В последующих выкладках мы будем использовать также диэлектрическую постоянную е =

Необходимо отметить, что к настоящему моменту существует немало работ, посвящённых поиску приближённых решений задач для уравнений гидродинамической модели (11), (12). Так, численное моделирование процесса переноса заряда для Ш задачи о баллистическом диоде рассматривается в [26]—[32], а поиск численных решений 2Э задач - в [33, 84]. Авторы работы [84] для создания непротиворечивой энергетически-транспортной модели переноса электронов в полупроводниковых устройствах используют смешанную схему конечных элементов.

Как мы видим, гидродинамическая МБР модель содержит квазилинейные нестационарные уравнения в частных производных. Таким образом, помимо проблемы конструирования вычислительных алгоритмов для поиска приближённых решений задач физики полупроводников, остро встаёт вопрос об обосновании построенных алгоритмов. Провести доказательство существования, единственности решения квазилинейных систем в общем случае не представляется возможным. В связи с этим, при создании вычислительных моделей необходимо чётко представлять, насколько найденные с помощью этих моделей приближённые решения будут соответствовать свойствам исходной дифференциальной задачи. На наш взгляд, правильное применение методов математического моделирования невозможно без детального анализа соотношения между исходной математической и вычислительной моделями явления. Заметим, что в теории дифференциальных уравнений уже давно существует метод одновременного изучения дифференциальной задачи и её конечно-разностного аналога (см. по этому поводу работу [35] и монографии [36, 38]).

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

Используя принципы, изложенные в [37], при конструировании и исследовании вычислительных моделей в данной диссертации мы будем отталкиваться от требования адекватности вычислительной модели исходной дифференциальной задаче. Под адекватностью будем понимать следующее: разностная модель строится так, чтобы с её помощью можно было доказать теорему существования решения исходной задачи. Подобные доказательства опираются на наличие разностных аналогов априорных оценок, которые удалось получить для уравнений МЕР модели в некоторых случаях (при дополнительных предположениях). Последнее обстоятельство представляется чрезвычайно важным фактом, поскольку благодаря наличию указанных оценок при проведении вычислений мы можем быть уверены в том, что приближённое решение в пределе действительно удовлетворяет свойствам дифференциальной задачи. Описанную методику, как мы увидим, можно реализовать при использовании различных численных методов.

К настоящему моменту создано впечатляющее многообразие способов построения вычислительных моделей для поиска решения различных задач. К ним относится и целый комплекс конечно-разностных методов построения явных, неявных схем, схем с весами, схем бегущего счёта (см. [39, 41, 42]); и метод прямых (см., например, [39, 51]); и методы частиц-в-ячейках (см. [43]); и вариационный подход, который приводит нас к задаче о минимизации функционалов в соответствующих пространствах (см. работы [44, 45]); и методы конечных элементов (см. [46]); и методы Монте-Карло для решения многомерных задач (см. [47]); и т.д.

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

Именно такие сложности возникают при поиске численных решений задач физики полупроводников, использующих уравнения МЕР модели. Проблема заключается в том, что функция плотности легирования р(х, у) и диэлектрическая постоянная е (см. (14)) принимают значения на несколько порядков меньшее, чем величины остальных параметров задачи, а значения напряжений смещения (их мы определим ниже при постановке конкретных задач) могут быть достаточно велики. Некоторые пути решения подобных проблем предложены, например, в книге А.Ф. Воеводина, С.М. Шугрина (см. [41. Ч. III, гл. 3]) и в работах А.И. Задорина при исследовании задач с погрничиым слоем (см. [48, 49]).

В данной диссертации мы предлагаем иной способ решения проблемы малых параметров - будем организовать дискретизацию таким образом, чтобы происходила автоматическая адаптация позиций узлов сетки к свойствам решения, погрешность аппроксимации была минимальной, а алгоритм, её использующий, давал высокую точность при небольшом количестве узлов. Подобная аппроксимация легла в основу вычислительного алгоритма, предложенного для поиска решений уравнений МЕР модели (13), (14) в гл. 1 данной диссертации. Вслед за автором описанной идеи К. И. Бабенко будем именовать схемы, автоматически учитывающие априорные свойства решения, схемами без насыщения или ненасыщаемыми схемами.

В [40] предложен весьма изящный способ построения вычислительных схем без насыщения, использующий для приближения неизвестных функций интерполяционный многочлен с узлами интерполяции в нулях полинома Чебышева. Подобные методы были впоследствии развиты в работах С.Д. Алгазина (см., например: [50]).

Но, тем не менее, описанный способ также не свободен от недостатков. Дело в том, что при использовании интерполяционных многочленов для аппроксимации неизвестной функции по каждой переменной мы в общем случае придём к линейным алгебраическим соотношениям с полностью заполненными матрицами. Работа с такими матрицами, как известно, влечёт за собой значительное увеличение временных затрат и рост погрешностей. Поэтому в диссертации для поиска численных решений уравнений МЕР модели в совокупности с интерполяционными многочленами мы будем использовать аппроксимацию сплайнами, которая тоже обладает сглаживающим эффектом и не допускает роста погрешностей (см. [71]). Кроме того, И.А. Задориным показано (см., например, [49]), что сплайн-интерполяцию можно модифицировать для расчёта задач с большими градиентами.

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

Основной целью данной диссертащш является конструирование и обоснование вычислительных моделей для поиска стационарных решений уравнений гидродинамической МЕР модели в 1Б и 2Б случаях. Кроме того, важнейшими задачами для нас являются создание на базе построенных моделей работоспособных вычислительных алгоритмов и их реализация на ЭВМ, а также апробация сконструированых вычислительных моделей при расчётах различных задач (в том числе проблем, не имеющих отношения к физике полупроводников) и сравнение найденных приближённых решений этих задач с результатами, полученными другими методами. Предлагаемые в данной работе алгоритмы основаны на применении метода установления.

Основная идея метода установления (см. [40]) состоит в том, чтобы искать стационарное (либо периодическое) решение нестационарной задачи как предел её решений при £ —> оо.

Замечание 2. При реализации подобной идеи мы будем использовать метод прямых и проводить дискретизацию уравнений МЕР модели по временной переменной, а затем осуществлять переходы с предыдущего временного слоя на следующий до тех пор, пока решение не «установится». Иначе говоря, алгоритм, использующий метод установления, будет остановлен только тогда, когда норма разности решений на следующем и предыдущем слоях станет достаточно близкой к нулю.

Прежде чем приступить к воплощению данной идеи, уместно задаться вопросом «В каких случаях мы можем гарантировать сходимость метода установления»? В книге [40] показано, что для сходимости метода достаточно, чтобы стационарное решение было асимптотически устойчивым по Ляпунову. В таком случае стремление к пределу будет происходить с экспоненциальной скоростью. Однако на практике добиться сходимости метода установления в расчётах с малыми параметрами (даже при наличии теоретических оценок, доказывающих эту сходимость) - задача нетривиальная. Для того чтобы справиться с ней, при проведении расчётов мы задействовали несколько методов: итерации по нелинейности, фильтр нелинейного сглаживания, применение прогонок в разных направлениях, «вытягивание» параметров и пр.

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

Подход 1

Стационарные уравнения МЕР модели

Ре1уляршашш

Регуляризованная нестационарная система

1 Метод прямых

Система ОДУ второго порядка и

Интерполяция сплайнами

Трёхточечиая разностная схема

Ц Метод прогонки

Решение регуляриэованных уравнений на каждом временном слое

Исходные нестационарные уравнения ¡13), {14)

Симметризация хг

Симметрическая по Фридрихсу система

Дискретизация ]

Класс "устойчивых" разностных схем

Преобразования хг

Трехточечнэя матричная схема

Метод матричной прогонки Д

Решение на каждом временном слое

Метод установления тг:

Стационарные решения

Рис. 1. Два подхода для поиска стационарных решений уравнений МЕР модели

Первая глава диссертации посвящена описанию подхода 1. Он базируется на сведении уравнений МЕР модели в стационарном случае к системе из трёх уравнений Пуассона (см. §1 гл. 1). В §2 гл. 1 мы рассмотрим модельную линейную краевую задачу для уравнения Пуассона

Дф = Я®, у) (15) здесь Ф(х,у) - неизвестная функция, f(x,y) - достаточно гладкая правая часть), на примере которой будем строить и обосновывать вычислительный алгоритм. Для (15) в §3 гл. 1 будет предложено три вида нестационарных регуляризаций: первый вид - параболическая регуляризация щ = Ди- f(x,у), второй вид - регуляризация Соболева ut - Aut = Аи- f(x, у), третий вид - гиперболическая регуляризация utt + Кщ - Дп + f(x, у) = 0.

Здесь u(t,x,y) - новая неизвестная функция, К > 1 - постоянная. Целесообразность подобного перехода от системы уравнений Пуассона к нестационарным регуляризованным уравнениям оправдана, потому что это позволило нам реализовать идею о поиске стационарного решения уравнений МЕР модели с помощью метода установления.

Обоснование метода установления в случае применения каждой из трёх регуляризаций представляет отдельную задачу, так как прежде всего необходимо доказать однозначную разрешимость регуляризованного уравнения, а затем асимптотическую устойчивость по Ляпунову его решения при t —> оо. В п. 3.1-3.3 §3 для каждого из трёх видов регуляризаций будут получены априорные оценки на нормы решений регуляризованных задач, из которых следует однозначная разрешимость и устойчивость стационарного решения.

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

При реализации подхода 1 в §4 гл. 1 мы используем метод прямых (см. [51, 52, 53, 54]) и принципы схем без насыщения. При этом аппроксимируем в регуляризованных уравнениях производную щ разностями, а вторую производную ихх приблизим с помощью интерполяционного многочлена с узлами интерполяции в нулях полинома Чебышева. В результате регуляризованные уравнения сводятся к системам ОДУ второго порядка. Важно отметить, что в случае использования параболической регуляризации для доказательства устойчивости полученной дифференциально-разностной модели из ОДУ нам удалось вывести разностный аналог априорной оценки, имеющей место для регуляризованных уравнений. Таким образом, по крайней мере в случае модельной линейной задачи для уравнения (15) мы обосновали адекватность вычислительной модели первого подхода.

Для поиска приближённых решений полученных систем ОДУ в §5 гл. 1 мы воспользуемся сплайн-интерполяцией. В результате придём к трёхточечной схеме, решение которой будем искать хорошо известным методом прогонки. Вычислительная схема подхода 1 апробирована при поиске численных решений тестовой краевой задачи для уравнения Пуассона с известным точным решением. В §6 гл. 1 исследована зависимость относительной погрешности численного решения тестовой задачи от различных параметров. Основные результаты гл. 1 опубликованы в работе [70].

Структура, заложенная в основу второго подхода, разработанного в гл. 2 для поиска стационарных решений уравнений МЕР модели (см. рис. 1), качественно отличается от идей подхода 1. Подход 2 базируется на технике построения разностных схем для гиперболических уравнений, предложенной в работах A.M. Блохина, Р.Д. Алаева. В монографии [37] рассматриваются смешанные задачи для симметрических по Фридрихсу t-гипсрболи-ческих систем с диссипативными граничными условиями (данные системы достаточно подробно изучены С. К. Годуновым в [36]).

Одна из рассмотренных в [37] задач следующая. Найти решение U симметрической ¿-гиперболической системы

AXJt + BUX -I- CUу = 0 в удовлетворяющее граничным условиям при х = 0:

U7 = SU11

17)

16) и начальным данным при t = 0: и(0,ж,у) = и0(ж,2/). 15

Здесь А > О, В, С - симметрические матрицы размера N х N, имеющие специфический вид (см. [37]): Щ\ ( Щ \ ( uNo+1 \

V = \J{t^y) = . , U = - , V"= - , UN J \ uNo J \ uNo+Nl J

S - вещественная постоянная матрица размерности vVq х , < < N -положительные целые числа, == {(i, x, у); t, x > 0, у G M1}.

Пусть граничные условия (17) являются строго диссипативными (см. [36, §15]):

4SU,U)|i=0>fco(U",U'OU, (19) где fco > 0 - постоянная. Тогда, следуя [36, §16], можно записать на гладких решениях системы (16) тождество интеграла энергии (при выполнении (19) это тождество называется диссипативным интегралом энергии):

AU, U) + ^(BU, U) + |-(CU, U) = 0, (20) с помощью которого в случае диссипативных граничных условий (17) для

16) может быть получена следующая априорная оценка:

J(t) < J(0), t > 0, (21)

00 где J{t) = f (f (^4U> U)dx)dy - интегральная норма решения U.

R 0

В работе [37] предлагается строить разностную модель для решения (16),

17), (18) так, чтобы она допускала наличие разностного аналога диссипа-тивного интеграла энергии (20). Это даст нам возможность получить энергетическую оценку (разностный аналог априорной оценки (21)), из которой будут следовать устойчивость предлагаемой схемы и адекватность сконструированной вычислительной модели исходной дифференциальной задаче.

Для реализации этих идей в рамках подхода 2 в §1 гл. 2 нестационарные квазилинейные уравнения МЕР модели (13) посредством введения симмет-ризатора будут преобразованы к симметрической по Фридрихсу ¿-гиперболической системе, схожей с (16), для которой в §2 будут получены тождество интеграла энергии и априорная оценка, аналогичные (20), (21) соответственно. Важно заметить, что эти тождество и оценка были получены благодаря наличию двух форм записи симметрической системы. Учитывая принципы монографии [37], в §3 мы будем строить разностные схемы для симметрической системы уравнений МЕР модели так, чтобы для них был выполнен аналог априорной оценки, доказанной для исходной дифференциальной задачи. Основные результаты гл. 2 опубликованы в [55], [56].

Идеи и методы первого и второго подходов использованы в данной работе для поиска приближённых решений нескольких прикладных задач. В гл. 3 мы апробируем разработанные подходы при построении алгоритмов поиска стационарных решений ID задачи о баллистическом диоде. В §1 гл. 3 будет дана постановка проблемы и показано, что она поставлена верно. Все нюансы использования подхода 1 для поиска решений задачи о баллистическом диоде приведены в п. 2.1-2.4 §2. В §3 для указанной задачи будут созданы вычислительные модели, основанных на методе ортогональной прогонки, технике интегральных уравнений и схеме предиктор-корректор. Параграф 4 гл. 3 по-свящён конструированию класса «устойчивых» разностных схем подхода 2 в 1D случае (см. п. 4.1) и созданию на основе полученных соотношений вычислительного алгоритма для поиска решений задачи о баллистическом диоде (см. п. 4.2). Далее, в §5 гл. 3 мы предложим технику преобразования стационарных уравнений 1D МЕР модели в случае е = 0 (е - диэлектрическая постоянная) к системе алгебраических соотношений.

Существенным достижением является то, что все сконструированные в гл. 3 алгоритмы удалось реализовать на ЭВМ с помощью языков программирования Delphi 6 (среда Object Pascal) и Java (все тонкости этого процесса описаны в п. 6.1-6.3 §6). Вычислительная схема подхода 2 была распараллелена для проведения быстрых расчётов на кластере НГУ (см. https://www.nusc.ru/) Схема распараллеливания приведена в п. 6.4. Другой важный итог проделанной в гл. 3 работы - сравнение эффективности разработанных подходов и других предложенных методов при различных значениях параметров задачи (см. п. 6.5 §6). Нужно отметить, что найти стационарное решение задачи о баллистическом диоде для некоторых величин параметров удалось только при использовании подходов 1 и 2. Подход 2 позволил получить хороший результат при варьировании значений параметров в широком диапазоне. Графики численных решений, полученные в ходе расчётов, представлены в п. 6.5. Основные результаты гл. 3 опубликованы в работах [58, 59].

Алгоритм подхода 1 хорошо зарекомендовал себя при поиске решения 1D задачи о баллистическом диоде (см. гл. 3) и 2D задачи о переносе заряда в транзисторе MESFET (см. по этому поводу [70]). Как мы увидим, предложенные в подходе 1 методы гармонично сочетаются и в совокупности представляют собой универсальный и удобный «комбайн», учитывающий идеи схем без насыщения и решающий проблему малых параметров, входящих в МЕР модель. В гл. 4 диссертации и приложениях мы применим этот «комбайн» при поиске решений прикладных задач. В связи с этим введём понятие технологии построения вычислительных моделей для поиска решений 1Б и 2Б нестационарных смешанных краевых задач (далее просто «технология»). В основу данного понятия мы положим следующую схему.

В Ш случае: аппроксимация производной неизвестной функции по времени разностным отношением —» сведение задачи к поиску решения ОДУ второго порядка с краевыми условиями —» использование сплайн-функций для интерполяции решения полученного ОДУ —» сведение задачи к трёхто-чечной схеме с граничными соотношениями —» поиск решения на каждом шаге по времени методом прогонки.

В 2Б случае: аппроксимация производной неизвестной функции по времени разностным отношением, а производной по одному из пространственных направлений - интерполяционным многочленом с узлами интерполяции в нулях полинома Чебышева —» сведение проблемы к краевой задаче для системы ОДУ второго порядка —> поиск решения полученной системы в виде интерполяционного кубического сплайна класса С2 —> вывод трёхточечных матричных соотношений с граничными условиями —» поиск решения на каждом временном слое методом матричной прогонки.

В гл. 4 мы применим разработанную технологию при поиске стационарных решений задачи о переносе заряда в транзисторе МОБРЕТ с напоканалом из оксида кремния (представление такого транзистора дано в [84]). Постановка данной задачи для уравнений МЕР модели (13), (14) приведена в §1. Здесь же, отталкиваясь от результатов гл. 1, мы записали соотношения (13), (14) в виде системы из трёх уравнений Пуассона, В §2 гл. 4 для электрического потенциала 99 на множестве, где наноканал примыкает к остальной части транзистора МОБРЕТ было выведено дополнительное краевое условие, что позволило нам применить для поиска решения поставленной задачи разработанную технологию. При использовании технологии в совокупности с регуляризацией методом установления и рядом вспомогательных средств, обеспечивших сходимость метода установления, было получено искомое стационарное решение. Схема построения вычислительной модели, подробно описанная в гл. 4, изображена на рис. 2.

Кроме того, для поиска стационарных решений задачи о переносе заряда в

-ГТ

Регулпри «щия (параболическая или СооолеБскип)

Нестационарные уравнения тг

Технология

Решение регуляризованных уравнений на п-м врел!енном слое

Метод уст,1 лиоплепнл

Стационарное решение задачи

Рис. 2. Схема вычислительной модели гл. 4 транзисторе МОБРЕТ мы использовали хорошо известный метод продольно-поперечной прогонки (п.п.п.). Описание вычислительной модели, созданной на его основе, представлено в §3 гл. 4. Сравнение эффективности и быстродействия алгоритмов, построенных с помощью технологии и метода п.п.п., приведено в §4 гл. 4. Также здесь размещены графики стационарных решений, нолученные при реализации указанных методов на ЭВМ. Основные результаты гл. 4 опубликованы в [57]

Предложенная нами технология апробирована в данной работе при поиске численных решений задач, не имеющих отношения к физике полупроводников. Так, в Приложении III изучен вопрос об устойчивости ударных волн в сжимаемом вязком газе, движение которого описывается известными уравнениями Навье-Стокса. Как и в работах [60, 61] рассмотрена следующая ситуация: равномерный сверхзвуковой набегающий поток отделяется от возмущённого течения вязкого газа ударной волной с уравнением х = /(¿, у), где £ - время, (х, у) - декартова система координат (в Приложении III мы ограничимся рассмотрением плоского случая).

В §1 Приложения III после перехода к криволинейной ортогональной системе координат а = х—^-, (3 = уех уравнения Навье-Стокса для сжимаемого теплопроводного иолитропного газа преобразованы к четырём соотношениям для поиска удельного объёма V и температуры газа Г, а также компонент вектора скорости газа уа, у р. Для этих четырёх уравнений из обобщённых соотношений Рэнкина-Гюгонио, записанных на фронте ударной волны, получены граничные условия. При исследовании вопроса об устойчивости ударных волн в Приложении III последовательно решаются две задачи.

1. В §2 рассматривается задача о поиске стационарных значений компонент скорости уа, Ур в окрестности линии (3 = 0 {у = 0), которая сводится к поиску решений ОДУ с условиями на границах а = 0 и а = оо. После применения к ОДУ параболической регуляризации, метода установления и разработанной технологии были получены искомые стационарные значения уа, Ур при достаточно «напряжённых» значениях параметров (схема конструирования вычислительной модели, использованная при этом, приведена на рис. 3,а).

2. Затем в §3 Приложения III осуществляется линеаризация четырёх нестационарных уравнений для поиска V, Т, уа, ур, которые рассматриваются для малых возмущений искомых величин. Найденные в окрестности линии /3 = 0 компоненты вектора скорости уа, ур входят в данную систему в качестве коэффициентов. Чтобы задействовать созданную технологию, прежде всего мы применили к полученным в результате линеаризации 2Б соотношениям преобразование Фурье. В итоге пришли к Ш нестационарной краевой задаче в области {£ > 0, о; > 0} с независимым параметром Фурье. Для нахождения стационарного решения этой задачи мы используем метод установления в совокупности с технологией. Схема поиска решения в данном случае изображена на рис. 3,6. а б

Стационарные уравнения для поиска оа, в окрестности прямой /? = 0

Параболическая• регуляризация

Нестационарные уравнения Л ехнологня

Найдено решение регуляризованных уравнений на каждом временном слое

Метод установления

Стационарное решение задачи

Нестационарные уравнения для поиска малых возмущений неизвестных

-Технология

Решение на каждом временном слое

МетодЦ установления

Стационарное решение задачи

Рис. 3. Схемы вычислительных моделей Приложения III: а - модель для поиска значений пп, Ур в окрестности [3 = 0; б - модель для поиска значений малых возмущений V, Т, уа, ир

В ходе численных экспериментов при определённых значениях параметров удалось получить нулевое стационарное решение описанной задачи для малых возмущений V, Т, уа, уд. Данное обстоятельство является косвенным доказательством устойчивости ударных волн в сжимаемом вязком газе при рассмотренных значениях параметров.

В Приложении IV рассмотрена другая задача, напрямую связанная с добычей нефти. Известно, что проблема повышения отдачи нефтяных пластов имеет важное значение для современной энергетики. Трудность решения этой проблемы заключается в том, что в процессе эксплуатации в трещиноватых зонах коллекторов формируются водонефтяные слоистые системы, которые, блокируя транспортную структуру коллекторов, выводят значительные нефтеносные области из режимов водного вытеснения. Восстановление проницаемости коллектора возможно лишь в условиях разрушения слоистых водонефтяных структур. В [62] (см. также обзор [63]) была предложена и подробно описана новая гидродинамическая модель для газосодержащих нефтяных слоистых систем. В [63] выявлены факторы, позволяющие управлять устойчивостью таких слоистых структур. При наличии начального градиента концентрации газовой фазы и периодических внешних акустических возмущений в распределённой слоистой системе возникает параметрический резонанс, приводящий к выходу газа и разрушению всей структуры (см. по этому поводу также работы [64]-[66]).

В работе [66] для численного исследования параметрического резонанса использовалась широко известная схема Годунова. Однако из-за проблем, связанных с аппроксимацией граничных условий, некоторые аспекты численных результатов из [66] подверглись авторами позже сомнению. В Приложении IV с использованием созданной технологии строится численный алгоритм, свободный от этих недостатков.

В §1 Приложения IV кратко описана система гидродинамических уравнений для газосодержащих нефтяных слоистых систем. С учётом весьма правдоподобных упрощающих предположений после преобразований данной системы для малых деформаций и слоистой структуры поставлена двумерная нестационарная смешанная задача. Для поиска численного решения полученной задачи в §2 использована технология (схема поиска решения в дан-пом случае изображена на рис. 4). В результате при определённых значениях параметров обнаружен эффект роста со временем амплитуды колебания решения. Это и есть искомый параметрический резонанс, приводящий к разрушению слоистой структуры.

Рис. 4. Схема вычислительной модели Приложения IV

Таким образом, на примере описанных задач мы убедились, что разработанная технология является эффективным, удобным, надёжным и универсальным инструментом поиска приближённых решений Ши2Б задач определённого класса. В совокупности с методом установления и применением нестационарных регуляризаций технология позволяет искать решения стационарных уравнений (см. рис. 2, 3,а); стационарные решения нестационарных уравнений (см. рис. 3}б) и исследовать динамику решений нестационарных систем (см. рис. 4). Кроме того, рассматриваемая технология основана на идеях схем без насыщения и принципе адекватности, что является её несомненным преимуществом.

На защиту выносятся следующие основные положения.

1. Разработаны и теоретически обоснованы два новых подхода для поиска стационарных численных решений уравнений гидродинамической МЕР модели:

- первый подход является оригинальным сочетанием набора регуляризаций, идей схем без насыщения, сплайн-интерполяций и метода установления;

- второй подход представляет новый способ конструирования и исследования устойчивых разностных схем для нелинейных задач.

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

3. При использовании технологии подхода 1 разработан вычислительный алгоритм для иоиска стационарных решений задачи о переносе заряда в транзисторе МОЭРЕТ.

4. Все алгоритмы, описанные в работе, реализованы на ЭВМ, проведены численные эксперименты, результаты которых представлены в диссертации.

Автор выражает глубокую благодарность и признательность своему научному руководителю д-ру физ.-мат. наук, профессору Александру Михайловичу Блохину за постоянное внимание и помощь в работе, канд. физ.-мат. наук, Алесе Сергеевне Ибрагимовой за консультации при реализации алгоритмов и ценные советы в ходе исследований, а также Станиславу Андреевичу Боярскому и Роману Евгеньевичу Семенко за обсуждения и поддержку при выполнении работы.

Итак, основные результаты диссертации состоят в следующем:

1. Для поиска стационарных численных решений уравнений гидродинамической МЕР модели разработаны и обоснованы два новых эффективных подхода.

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

3. На основе идей подхода 1 создана универсальная технология конструирования вычислительных моделей для поиска решений 1D и 2D нестационарных смешанных краевых задач.

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

5. Все алгоритмы, описанные в работе, реализованы на ЭВМ с помощью языков Delphi б (среда Object Pascal) и Java.

Данная диссертация представляет собой исчерпывающее исследование поставленных задач. Однако в ней также обозначен круг проблем, требующих, по мнению автора, дальнейшего изучения. Так, класс «устойчивых» разностных схем, предложенный в гл. 2 работы, может быть использован для поиска приближённых решений 2D задач физики полупроводников, например задач о переносе заряда в транзисторах MESFET (подробное описание данного прибора приведено в [25], [33]) и MOSFET. Для строгого обоснования устойчивости разностных схем необходимо доказать, что сделанные в §2, 3 гл. 2 предположения относительно ограниченности норм матриц и функций, входящих в симметрическую систему, а также разностные аналоги этих предположений, действительно выполняются. Технику преобразования уравнений 1D МЕР модели при е = 0, подробно описанную в §5 гл. 3, можно перенести на случай, когда £ ^ 0. Как было показано в работах Э.А. Бибердорф, Н.И. Поповой (см., например, [73]), метод ортогональной прогонки позволяет проводить вычисления с гарантированной точностью. Организация подобных вычислений и аккуратная оценка точности при поиске решений задачи о баллистическом диоде являются актуальными задачами.

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

Заключение

В гл. 1, 2 данной работы, используя метод установления, мы разработали два подхода для поиска стационарных решений уравнений гидродинамической МЕР модели. Первый подход основан на сведении уравнений модели к регуляризованной системе. При поиске её приближённых решений была создана эффективная и удобная технология, основанная на идее метода прямых, схем без насыщения и принципе адекватности. Второй предложенный подход олицетворяет новую технику конструирования и анализа «устойчивых» разностных схем на нелинейном уровне. Данный подход опирается на наличие двух симметрических форм записи исходной системы, принцип адекватности и ряд дополнительных предположений.

Высокая эффективность описанных подходов продемонстрирована в гл. 3 при поиске стационарных решений задачи о баллистическом диоде. В последующих разделах диссертации (а именно, в гл. 4 и Приложениях III, IV) сделан акцент на технологии конструирования вычислительных моделей для поиска решений и 2D нестационарных смешанных краевых задач, предложенной в гл. 1. При использовании данной технологии в гл. 4 найдены стационарные решения задачи о переносе заряда в транзисторе МОБЕЕТ, а в Приложениях III, IV технология апробирована при исследовании прикладных проблем, не имеющих отношения к физике полупроводников.

Важно отметить, что область применимости созданной технологии достаточно широка - она ограничена требованием наличия априориой гладкости разыскиваемого решения, размерностью задачи (Ю либо 2Б) и конфигурацией множества, в котором поставлена задача (для 2Б уравнений - прямоугольное множество, для 1Б - прямая). Однако в случае применения технологии в совокупности с методом установления при поиске стационарных решений прикладных задач возникают трудности с обеспечением сходимости данного метода. В диссертации предложено несколько способов решения этой проблемы, но, тем не менее, в некоторых случаях добиться сходимости не удаётся.

1. Roosbroek W. van. Theory of flow of electrons and holes in germanium and other semiconductors // Bell System Techn. J. 1950. Vol. 29. P. 560-607.

2. BEIRAO DA VEIGA H. On the semiconductor drift-diffusion equations // Diff. Integral Equat. 1996. Vol. 9. No. 4. P. 729-744.

3. Gajewski H., groger K. Semiconductor equations for variable mobilities based on Boltzman statistics for Fermi-Dirac statistics // Math. Nachr. 1989. Vol. 140. P. 7-36.

4. Jin Liang. On a nonlinear integrodifferential drift-diffusion semiconductor model // SIAM J. Math. Anal. 1994. Vol. 25. No. 5. P. 1375-1392.

5. MOCK M.S. An initial value problem from semiconductor device theory / / SIAM J. Math. Anal. 1974. Vol. 5. No. 4. P. 597-612.

6. Mock M.S. On equations 'describing steady-state carrierdistributions in semiconductor device // Commun. Pure Appl. Math. 1972. Vol. 25. No. 6. P. 781-792.

7. Selberherr S. Analysis and Simulation of Semiconductor Devices. Wien; N. Y.: SpringerVerlag, 1984.

8. Grasser Т., Tang T.-W., Kosina H., Selberher S. A rewiew of hydrodynamic and energy-transport models for semiconductor device simulation // Proc. IEEE 91. 2003. No. 2. P. 251-274.

9. HAILIANG Li, MARKOWICH P. A. A review of hydro dynamical models for semiconductors: asymptotic behavior // Boletim da Sociedade Brasileira de Mathematica. 2001. Vol. 32. No. 3. P. 321-342.

10. AbdallaH N. В., Degond P. On a hierarchy of macroscopic models for semiconductors // J. Math. Phys. 1996. Vol. 37. No. 2. P. 3308-3333.

11. BLOTEKJAER K. Transport equations for electrons in two-valley semiconductors // IEEE Trans. Electron Devices. 1970. Vol. ED-17. P. 38-47.

12. Baccarani G., Wordman M. R. An investigation on steady-state velocity overshoot in silicon // Solid-State Electronics. 1982. Vol. 29. P. 970-977.

13. Li H., Markowich P., Mei M. Asymptotic behaviour of solutions of the hydrodynamical model of semiconductors // Proc. of the Royal Society of Edinburg. 2002. Vol. 132A.1. P. 359-378.

14. Gamba I. M. Stationary transonic solutions of a one-dimensional hydrodynamic model for semiconductors // Commun. PDE. 1992. Vol. 17. No. 3-4. P. 553-577.

15. Gardner C. L., Jerome J. W., Rose D. J. Numerical methods for the hydrodynamic device model: subsonic flow // IEEE Transactions on Computeraided Design. 1989. Vol. 8. P. 501-507.

16. Gardner C.L. Numerical simulation of a steady-state electron shock wave in asubmicrometer semiconductor device // IEEE Transactions on Electron Devices. 1991. Vol. 38. No. 2. P. 392-398.

17. HsiAO L. Some mathematical analysis in semiconductor devices // Methods and Applications of Analysis. 2001. Vol. 8. No. 4. P. 635-644.

18. ROMANO V. Nonparabolic band hydrodynamical model of silicon semiconductors and simulation of electron devices // Math. Meth. Appl. Sci. 2001. Vol. 24. P. 439-471.

19. Anile A. M., Junk M., Romano V. et al. Cross-validation of numerical schemes for extended hydrodynamical models of semiconductors // Math. Models Meth. Appl. Sci. 2000. Vol. 10. P. 833-861. .

20. Anile A. M., Romano V., Russo G. Extended hydrodynamical model of carrier transport in semiconductors // SIAM J. Appl. Math. 2000. Vol. 61. P. 74-101.

21. Blokhin A. M., Iordanidy A. A. Numerical investigation of a gas dynamical model for charge transport in semiconductors // COMPEL. 1999. Vol. 18. P. 6-37.

22. Romano V., Russo G. Numerical solution for hydrodymamical models of semiconductors // Math. Models Meth. Appl. Sei. 2000. Vol. 10. P. 1099-1120.

23. Anile A. M., Romano V. Hydrodynamical modeling of charge carrier transport in in semiconductors // Meccanica. 2000. Vol. 35. P. 249-296.

24. Romano V. 2D simulation of a silicon MESFET with a non-parabolic hydrodynamical model based on the maximum entropy principle //J. Comp. Phys. 2002. Vol. 176. P. 70-92.

25. Romano V. 2D Numerical Simulation of the МЕР Energy-Transport Model with a Finite Difference Scheme // J. Comp. Phys. 2007. Vol. 221. P. 439-468.

26. Курант P., Фридрихс К., JIebh Г. О разностных уравнениях математической физики // Успехи мат. наук. 1940. Вып. 8. С. 125-160.

27. Годунов С. К. Уравнения математической физики. М.: Наука, 1979. 392 с.

28. БЛОХИП А. М., АЛАЕВ Р. Д. Интегралы энергии и их приложения к исследованию устойчивости разностных схем. Новосибирск, 1993.

29. Ладыженская О. А. Краевые задачи математической физики. М.: Наука, 1973.

30. Самарский А. А. Теория разностных схем: Учеб. пособие для вузов. 3-е изд., испр. М.: Наука, 1989. 616 с.

31. ЛАЕВСКИЙ Ю. М. Методы конечных элементов (основы теории, задачи). Новосибирск, 1999. 166 с.

32. Березин И. е., Жидков Н. п. Методы вычислений. М.: Физматгиз, 1962. Т. 2. 620 с.

33. Rothe Б. Н. Zweidimensionale parabolische Randwertaufgaben als Grenzfall eindimensionaler Randwertaufgaben // Math. Ann. 1929. Bd. 102, Heft 4/5.

34. Будак Б. M. К решению краевой задачи параболического типа // Вестн. МГУ. 1955. № 8. С. 33-38.

35. Крылов В. И., Лискавец О. А. Оценка погрешности метода прямых для задач Гурса // ДАН БССР. 1963. Т. 7. № 8. С. 505-509.

36. Blokhin А. М., Boyarsky S. A., Semisalov В. V. On an approach to the construction of difference schemes for the moment equations of charge transport in semiconductors // Le Matematiche. 2009. Vol. LXIV. Fase. I. P. 77-91.

37. Блохин A. M., Ибрагимова А. С., Семисалов Б. В. Конструирование вычислительных алгоритмов для задачи о баллистическом диоде // Выч. мат. и мат. физ. 2010. Т. 50. № 1. С. 188-208.

38. Блохин А. М., Семисалов Б. В. Конструирование одного класса вычислительных схем в задаче о баллистическом диоде // Мат. моделирование. 2010. Т. 22. № 7. С. 3-21.

39. Blokhin А. М. On stability of shock waves in a compressible viscous gas // Le Matematiche. 2002. Vol. LVII. Fase. I. P. 3-19.

40. Blokhin A. M., Trakhinin Yu. L. On a modified shock front problem for the compressible Navier-Stokes equations. // Quar. Appl. Math. 2004. Vol. 62. P. 222-234.

41. Dorovsky V. N., Dorovsky S. V. A hydrodynamic model of water-oil layered systems containing gas // Math. Comput. Mod. 2002. Vol. 35. P. 751-757.

42. Белоносов B.C., Доровский В.Н., Белоносов А.С. и др. Гидродинамика га-зосодержащих слоистых систем // Успехи механики. 2005. Т. 3. № 2. С. 37-70.

43. Dorovsky V. N., belonosov V. S., belonosov A. S. Numerical investigation of parametric resonance in water-oil structures containing gas // Math. Сотр. Mod. 2002. Vol. 36. P. 203-209.

44. Доровский С. В., Доровский В. Н. О возможностях методов электроразведки в исследовании устойчивости водонефтяных слоистых систем // Геология и геофизика. 2006. Т. 47. № 7. С. 892-901.

45. Блохип А. М., Доровский С. В., Доровский В. Н. О возможностях методов электроразведки в исследовании устойчивости водонефтяных слоистых систем. Часть II // Геология и геофизика. 2006. Т. 47. № 11. С. 1-21.

46. ГаевскиЙ X., Грегер К., Захариас К. Нелинейные операторные уравнения и операторные дифференциальные уравнения. М.: Мир, 1978.

47. Блохин А. М., Ибрагимова А. С., Красников Н. Ю. Об одном варианте метода прямых для уравнения Пуассона // Выч. технологии. 2007. Т. 12. № 2. С. 33-42.

48. Белман Р. Введение в теорию матриц. М.: Наука, 1976.

49. Блохин А. М., Ибрагимова А. С., Семисалов Б. В. Конструирование вычислительного алгоритма для системы моментных уравнений, описывающих перенос заряда в полупроводниках // Мат. моделирование. 2009. Т. 21. № 4. С. 15-34.

50. Завьялов Ю. С., Квасов Б. И., Мирошниченко В. JI. Методы сплайн-функций. М.: Наука, 1980.

51. Кузнецов С. В. Решение краевой задачи для обыкновенных дифференциальных уравнений // Тр. ИМ СО АН СССР. Выч. методы лин. алгебры. Новосибирск: Наука, 1985. Т. 6. С. 85-110.

52. Бибердорф Э. А., Попова Н. И. Контроль точности решения краевой задачи методом ортогональной прогонки. Новосибирск, 2009 (Препринт / РАН. Сиб. отд-ние. ИЯФ; № 1).

53. Блохин А. М., соковиков И. Г. Об одном подходе к конструированию разностных схем для квазилинейных уравнений газовой динамики // СМЖ. 1999. Vol. 40. No. 6. P. 1236-1243.

54. Blokhin A. M., Bushmanov R. S., Romano V. Nonlinear asymptotic stability of the equilibrium state for the МЕР model of charge transport in semiconductors // Nonlinear Analysis. 2006. Vol. 65. P. 2169-2191.

55. Blokhin A. M., Bushmanov R. S., Romano V. Global existence for the system of the macroscopic balance equations of charge transport in semiconductors //J. Math. Anal. Appl. 2005. Vol. 305. R 72-90.78