ИНСТИТУТ МАТЕМАТИЧЕСКОГО МОДЕЛИРОВАНИЯ
РОССИЙСКОЙ АКАДЕМИИ НАУК
На правах рукописи
Корякин
Павел Владимирович
Моделирование полупроводникового диода
Специальность 05.13.18 – математическое моделирование,
численные методы и комплексы программ
Автореферат
диссертации на соискание ученой степени
кандидата физико-математических наук
Москва
2010
Работа выполнена в Институте математического моделирования РАН
Научный руководитель – к.ф.-м.н. Е.А. Альшина
Официальные оппоненты:
– доктор физико-математических наук, профессор Боголюбов Александр Николаевич
– доктор физико-математических наук, профессор Тупчиев Виль Асадулаевич
Ведущая организация – Физический факультет Московского Государственного Университета.
Защита состоится «____» _______________________________2010 года в _______часов на заседании диссертационного совета Д 501.002.09 при Московском государственном университете имени М.В. Ломоносова по адресу: 119991, Москва, Ленинские горы, дом 1, стр. 4, НИВЦ МГУ.
С диссертацией можно ознакомиться в библиотеке ИММ РАН.
Автореферат разослан «_____»________________2010г.
Ученый секретарь диссертационного совета
кандидат физико-математических наук Змитренко Н.В.
Общая характеристика работы
Актуальность темы. Диффузионно-дрейфовая модель является самым популярным упрощением фундаментальной системы уравнений полупроводника [1-2]. Модель состоит из пяти уравнений, которые задают взаимосвязь между пятью основными характеристиками диода: концентрации электронов и дырок, дырочный и электронный токи, а также напряженность электрического поля.
Диффузионно-дрейфовая модель лежит в основе большинства современных моделей, отличия которых друг от друга и от родительской модели заключаются, как правило в том, что их разработчики, либо пытаясь учесть особенности конкретных приборов, либо преодолевая проблемы счёта, изменяют многочисленные коэффициенты, делают их нелинейными и так далее. Но все эти модификации не меняют базу модели.
Два из пяти уравнений модели являются дифференциальными уравнениями первого порядка по времени, и все пять уравнений являются дифференциальными уравнениями первого порядка по пространству. Таким образом, модель требует постановки пяти граничных условий. Изученные работы, в частности [1-2], посвященные моделированию диодов как правило предлагают задавать граничные условия следующим образом: четыре условия получаются фиксированием концентраций дырок и электронов на обоих контактах прибора, а в качестве пятого условия берётся естественное равенство разности потенциалов на границах и интеграла от напряженности поля. Реже предлагается полагать равной нулю напряженность поля на границах (из тех соображений, что вне полупроводника, в металле, падение напряжение очень мало, по сравнению с падением его в диоде), а также накладывается какое-либо условие на токи. Пятое, интегральное условие при этом сохраняется.
Однако, проведённый в работе анализ показал, что эти условия не только не вполне оправданы физически, но также при их использовании возникает проблема несовместимости граничных условий с самой системой. Так, например, возникает широко известная проблема, когда уравнение наряду с диссипативным членом содержит ещё и конвективный, что сильно усложняет задачу постановки условий, а точнее задачу численного интегрирования со вроде бы математически правильными условиями.
В силу целого ряда причин, численные расчёты по диффузионно-дрейфовой модели являются довольно сложной задачей. Как и большинство проблем микроэлектроники, рассматриваемую модель чаще всего приходится применять либо в слоистых средах, в которых свойства материалов меняются скачкообразно, либо в средах с непрерывно, но очень быстро меняющимися свойствами. Кроме того, сами величины, участвующие в модели могут иметь очень большие градиенты. Описанные выше проблемы приводят к тому, что систему уравнений модели очень трудно интегрировать численно, трудно построить устойчивый численный алгоритм. Практики как правило сводят систему к системе меньшей размерности, путём введения некоторых упрощений, что даёт им возможность интегрировать эту систему явными схемами. Кроме того, система уравнений рассматриваемой модели , как и многие другие задачи микроэлектроники имеет особенности, то есть точное решение может быть не гладким или даже не непрерывным, что сильно усложняет задачу численного интегрирования и требует применения особенно надёжных численных методов.
Цель работы - во-первых, разработать надёжный численный алгоритм, позволяющий решать широкий круг задач в слоистых средах, написать разностную схему для интегрирования системы уравнений диффузионно-дрейфовой модели. Во-вторых, разработать методику, позволяющую при численном интегрировании дифференциальных уравнений контролировать точность получаемого решения, диагностировать наличие и определять тип особенности точного решения. Разработать надёжный численный алгоритм, для интегрирования системы уравнений модели, с возможностью задания различных вариантов граничных условий.
Научная новизна. Впервые была предложена простая и очень эффективная методика диагностики особенностей точных решений при численном интегрировании ОДУ. Предложен новый тип разностных схем – так называемые бикомпактные схемы, применение которых особенно актуально для решения задач в слоистых средах. Пространственная аппроксимация бикомпактного типа записана для системы уравнений диффузионно-дрейфовой модели полупроводникового диода.
Практическая ценность работы. Разработан надежный численный алгоритм для интегрирования системы уравнений модели полупроводникового диода, реализованный на встроенном языке математического макета Matlab. Он позволяет задавать характеристики моделируемого прибора, ставить различные варианты граничных условий и проводить расчет статических и динамических характеристик диода. Кроме того, программа автоматически контролирует точность расчета по методу Ричардсона и анализирует наличие сингулярностей у точного решения. В случае обнаружения сингулярности, программа диагностирует её тип.
Апробация работы. Полученные результаты докладывались и обсуждались на нескольких российских и международных конференциях, среди которых были Международный конгресс математиков в Мадриде в 2006 года, конференция памяти А.Ф. Сидорова «Актуальные проблемы прикладной математики и механики», Всероссийская школа-семинар “Современные проблемы математического моделирования”. По материалам диссертации сделан доклад на совместном семинаре Института математического моделирования РАН и кафедры математического моделирования Московского физико-технического института (май 2008).
Структура и объем работы. Диссертация состоит из введения, трех глав, заключения, списка литературы и приложения. Общий объем диссертации 75 страниц, рисунков 34, таблиц 8. Список литературы включает 40 наименований.
Публикации. По теме диссертации всего опубликовано 14 работ, основные из которых представлены в конце автореферата.
Краткое содержание работы.
Введение включает описание диффузионно-дрейфовой модели полупроводникового диода. Описывается моделируемый прибор, рассматривается система уравнений модели и обсуждаются возникающие при численном интегрировании системы счётные проблемы. Также даётся небольшой научно-популярный экскурс в процесс производства полупроводников. Введение также содержит формулировку основных целей работы и краткое содержание глав.
Первая глава посвящена новому типу разностных схем – так называемым бикомпактным схемам. Сначала рассматриваются классические подходы к записи разностных аппроксимаций пространственных производных [3]. Обсуждаются проблемы, возникающие при использовании классических пространственных шаблонов, производится краткий анализ применимости подходов для задач в слоистых средах, каковыми являются большинство задач микроэлектроники. Далее обсуждаются компактные разностные схемы [4], приводится обзор результатов, полученных другими исследователями в этой области [5-7] и предлагается новый тип разностных схем, так называемые бикомпактные схемы.
Решение задач в слоистых средах осложняется тем, что очень трудно построить аппроксимацию, дающую высокий порядок точности на стыках сред. Если задавать сетку так, что граница сред лежит между узлами сетки, то построить адекватную аппроксимацию тяжело. Будем выбирать сетку так, чтобы её узлы попадали на границы слоёв. Такие сетки называют специальными. Однако если шаблон аппроксимации содержит три и более пространственных узла, то специальные сетки не спасают: когда внутренний узел шаблона совпадает с границей, приходится использовать аппроксимацию производных через разрыв коэффициентов. В этих узлах аппроксимация имеет пониженный порядок точности, что кардинально ухудшает общую точность расчёта. Использование в расчёте полуцелых узлов приводит к такому же эффекту.
Отметим здесь, что конкретно для уравнения теплопроводности А.А. Самарским была построена общеизвестная классическая схема, с аппроксимацией на трёхточечном шаблоне, то есть с использованием полуцелых узлов. Однако для уравнений других типов не получается построить хорошие схемы для слоистых сред с использованием полуцелых узлов. Оказалось, что для того, чтобы избежать описанных проблем, необходимо замкнуть всю схему в одном интервале и изгнать из расчетов полуцелые точки.
Бикомпактными называют схемы, шаблон которых состоит лишь из двух узлов сетки. Такой тип схем является мало исследованным, но при этом очень удобным в ряде практических задач. Такие схемы позволяют избавиться от некоторых нежелательных эффектов, как, например, эффект отражения уходящей волны от бесконечно удалённой жёсткой границы при расчетах схемами с использованием полуцелых узлов. При использовании бикомпактных схем этот нефизичный эффект исчезает. Расчёты с её использованием показали, что эффект отражения исчезает.
В качестве примера, иллюстрирующего методологию построения схем бикомпактного типа, приводится вывод двух бикомпактных схем для уравнения теплопроводности. Рассматривается следующая задача:
Она решается методом прямых: пространственные производные заменяются разностными аппроксимациями, затем полученная система дифференциальных уравнений интегрируется по времени.
Чтобы получить аппроксимацию второй производной точности  нужен 3-точечный шаблон. Однако наша цель построить бикомпактную схему, в которой шаблон по пространству будет состоять из двух точек, то есть задачу мы будем решать в пределах одного интервала сетки по пространству. Чтобы воспользоваться 2-точечным шаблоном, заменим эквивалентной системой двух уравнений первого порядка:
Помимо температуры  , здесь появляется тепловой поток  , который мы для удобства записи берём со знаком минус. Граничные условия задаются только для температуры - поток в граничных точках не известен. В такой постановке, помимо начального профиля температуры также необходимо задавать начальный профиль потока.
В неподвижной слоистой среде коэффициент  и свободный член  имеют неподвижные разрывы. Между разрывами считаем их многократно непрерывно дифференцируемыми. Благодаря наличию разрывов решение будет обобщённым. При этом физически правильным является решение, в котором  ,  всюду непрерывны.
Введём по пространству специальную сетку так, чтобы все точки, в которых функции и имеют разрывы, являлись бы узлами сетки. Схему для узловых значений построим методом прямых, интегрируя по пространству:
При этом под узловыми значениями функции надо подразумевать всегда односторонние пределы изнутри данного интервала. Но для теплопроводности величины  и  непрерывны всюду, в том числе на разрывах коэффициентов. Поэтому для них это просто узловые значения.
Задача получения бикомпактной пространственной аппроксимации заданной точности сводится к взятию интегралов в правых частях , с нужной точностью. Так, для получения схемы точности  , интегралы берутся по формуле трапеций, а для получения точности  , интегралы аппроксимируются по формуле Симпсона 4-го порядка.
В схеме второго порядка точности по пространству, алгебраическими преобразованиями удаётся исключить поток, в результате чего получается схема, которую формально можно назвать трёхточечной, однако она полностью эквивалентна двухточечной и сохраняет аппроксимацию для слоистых сред. Получаемая схема отличается от традиционной схемы для уравнения теплопроводности тем, что в левой части вместо производной по времени в центральном узле шаблона стоит линейная комбинация производных в трёх узлах. Это, незначительное на первый взгляд, отличие приводит к тому, что полученная схема обладает уникальным спектром -  [9]. Очевидно, что собственные значения этой схемы растут много быстрее, чем СЗ классической схемы -  и даже чем СЗ точного решения  .
Рис. 1. Спектры схем.
Рисунок 1 иллюстрирует поведение спектров рассматриваемых схем. Видно, что для бикомпактной схемы второго порядка высокие гармоники затухают гораздо быстрее не только классической схемы, но и точного решения, что обеспечивает получение более гладкого решения.
Первая глава заканчивается описанием тестовых расчётов, на которых исследовались свойства полученных схем и проводилось сравнение с классической схемой. Там же приводится анализ структуры ошибки численного решения и приводятся интересные рассуждения о природе нечетных степеней шага в разложении ошибки.
Во второй главе описывается методика диагностики особенностей точных решений при численном интегрировании дифференциальных уравнений. В общем виде, задачу, которую помогает решать предлагаемая методика можно сформулировать так. Предположим, что существует некоторый алгоритм решения определённого класса задач. Предположим так же, что есть конкретная задача, точное решение которой имеет особенность в какой-либо точке. Возникает вопрос: можно ли по результатам работы численного алгоритма понять, в какой точке решение имеет особенность и какого рода эта особенность?
Такая задача была поставлена в 2003 году на семинаре академика Г.И. Марчука в Институте вычислительной математики РАН. Тогда был предложен следующий пример: задача Коши для обыкновенного дифференциального уравнения (ОДУ):
 .
Точное решение задачи выглядит следующим образом:
 .
Оно имеет особенность в точке  и не существует при  .
Будем численно решать задачу (1.5), например, по явной схеме Эйлера порядка точности  с постоянным шагом 
 ,
где  - значение в следующий момент времени.
Численное решение положительно, монотонно возрастает и существует при сколь угодно больших  . По его виду не возможно сделать вывод о наличии полюса у точного решения. Кажется, что точное решение быстро возрастает и существует при любых  . Такое же качественное поведение дают явные схемы Рунге-Кутта более высоких порядков точности.
В данной работе приводится методика диагностики решения на наличие особенностей. Она основана на приеме сгущения сеток и позволяет выявлять не только сингулярности решения, но и более тонкие особенности: например, ограниченность числа непрерывных производных.
Методика основана на одностадийной схеме Розенброка с комплексным коэффициентом [8] и методе апостериорной оценки точности, предложенной Ричардсоном ещё в 1927 году.
Если проинтегрировать задачу (1.5) численно популярными явными схемами Рунге-Кутта различных порядков точности, чисто неявной схемой Розенброка и схемой Розенброка с комплексных управляющих коэффициентом, то результаты расчетов будут кардинально отличаться друг от друга на качественном уровне. Рисунок 2 иллюстрирует эти различия. Видно, что счет по схемам Рунге-Кутта (ERK1, ERK2, ERK4) разваливается сразу за моментом времени, в котором точное решение имеет особенность. Причем момент, в который происходит переполнение численного решения, не совпадает с моментом разрушения точного решения и невозможно сделать вывод даже о местоположении особенности. Решение по чисто неявной схеме Розенброка (обратный Эйлер) вблизи момента разрушения точного решения, пытаясь повторить гиперболу, уходит в область отрицательных значений, что не физично. И только решение по комплексной схеме Розенброка (CROS) ведёт себя уникальным образом – оно стабилизируется на некотором значении  сразу за моментом  . Причем высота этой полочки () зависит от числа узлов сетки, на которой ведется расчёт – чем подробнее сетка, тем выше полочка.
Рис. 2. Поведение численного решения для различных схем
Рассмотрим задачу (1.5) в более общем виде
 .
Её точное решение выглядит следующим образом -  . Применительно к (1.7) семейство одностадийных схем Розенброка выглядит следующим образом:
 .
При использовании схемы CROS с  принимает вид
 .
Из вида следует несколько легко доказуемых утверждений.
Утверждение 1. Существует значение , при котором численное решение схемы CROS не меняется при переходе на следующий временной слой  :  .
Утверждение 2 (свойство притяжения). При  приращение функции  и, наоборот, при  приращение функции  . Т.е. при любом значении  на текущем временном слое, следующий шаг схема CROS делает по направлению к положению равновесия.
Также удаётся описать условия устойчивости этого положения равновесия и сформулировать следующее утверждение:
Утверждение 3. При особенности точного решения ОДУ типа полюса  численное решение схемы CROS стремится к пределу  . При особенности типа корня  численное решение схемы CROS стабилизируется в полосе вокруг положения равновесия .
Для логарифмической особенности, когда точное решение ведет себя как  численное решение по схеме CROS также стабилизируется и аналогично случаю степенной особенности удаётся доказать следующие утверждения.
Утверждение 43. Существует положение равновесия  , при котором численное решение схемы CROS не меняется при переходе на следующий временной слой.
Утверждение 5. При приращение функции и, наоборот, при приращение функции .
Утверждение 6. В случае логарифмической особенности положение равновесия схемы CROS асимптотически устойчиво.
Таким образом, было установлено, что для рассмотренных выше типов особенностей численное решение стабилизируется на некотором значении .
Дельнейшие исследования показали, что учитывая тот факт, что значение зависит от величины шага сетки, на которой ведется расчёт и от характеристик точного решения, можно построить алгоритм, основанный на методе сгущения сеток Ричардсона, позволяющий диагностировать тип особенностей. Хорошо известно, что проведя расчёты на трёх вложенных сетках, сгущая их так, чтобы часть узлов более подробной сетки являлась бы узлами менее подробной, можно посчитать эффективный порядок точности численного решения. Если точное решение задачи обладает гладкостью, достаточной для реализации численным методом своего теоретического порядка точности, то эффективный порядок точности должен стремиться к теоретическому при увеличении числа узлов сетки. Оказалось, что можно доказать следующие два утверждения:
Утверждение 7. При наличии сингулярности точного решения типа  в узлах  эффективный порядок точности схемы CROS  при  .
Утверждение 8. При наличии сингулярности точного решения типа  в узлах погрешность схемы CROS  , а  при .
Основываясь на утверждениях 7 и 8 можно создавать программы, проводящие интегрирование ОДУ на сгущающихся сетках с автоматическим контролем точности и диагностикой особенностей. В самом простом варианте такая программа должна провести расчет на трех вложенных сетках и вычислить эффективный порядок в каждой точке первой сетки. Если начиная с какой-то точки эффективный порядок точности перестаёт быть равным теоретическому, то с точностью до шага исходной сетки можно сказать, что вблизи этой точки имеет место быть особенность точного решения. Более детальный анализ поведения эффективного порядка точности за обнаруженным моментом сингулярности позволит определить тип особенности.
Используя рекуррентное уточнение решения можно диагностировать более тонкие особенности решения, а именно, ограниченность числа непрерывных производных. Известно, что ошибку, вычисленную по методу Ричардсона по решениям на двух вложенных сетках можно учесть в качестве поправки для решения на исходной сетки, повысив тем самым порядок точности решения. Однако то, на сколько порядков удастся повысить порядок, зависит от свойств точного решения, а именно от его дифференцируемости. В работе приводится пример, иллюстрирующий одновременную диагностику как разрыва первой производной, так и следующей за ней степенной особенности.
В третьей главе выводится бикомпактная схема для системы уравнений диффузионно-дрейфовой модели. Как уже отмечалось выше, система уравнений модели состоит из пяти дифференциальных уравнений первого порядка по пространству, два из которых так же являются дифференциальными уравнениями по времени. Система интегрируется по методу прямых, то есть пространственные производные аппроксимируются разностными аналогами, в результате чего получается система обыкновенных дифференциальных уравнений по времени, которая интегрируется одностадийной схемой Розенброка с комплексным коэффициентом, на уникальные свойства которой мы неоднократно указывали.
Построение бикомпактной пространственной аппроксимации происходит следующим образом: вводится некоторая сетка по пространству и в каждом её интервале все уравнения системы интегрируются по  . Часть выражений, в которых стоят производные по удаётся проинтегрировать точно, остальные же интегралы аппроксимируются по формуле трапеций. В результате для каждого интервала сетки получается система уравнений, приведённая ниже:
Полученная система уравнений интегрируется по времени схемой Розенброка с комплексным коэффициентом CROS.
Программа, реализующая алгоритм интегрирования, описанный выше также может проводить расчет на сгущающихся сетках с автоматическим контролем точности по каждой переменной и диагностикой особенностей точного решения.
В заключении сформулированы основные результаты работы.
В приложении приводится текст программы расчета статических и динамических характеристик полупроводникового диода.
Основные результаты
Построен и исследован новый тип разностных схем применительно к уравнению теплопроводности. Построены схемы разных порядков точности. Исследована устойчивость схем. Исследована структура ошибки численного решения для разных видов сеток. Теоретические проработаны подходы построения бикомпактных схем для двумерных задач.
Разработана оригинальная методика диагностики особенностей точных решений при численном интегрировании обыкновенных дифференциальных уравнений, основанная на методе апостериорной оценки точности Ричардсона и применении одностадийной схемы Розенброка с комплексным коэффициентом.
Подробна описана методика написания программ интегрирования ОДУ с контролем точности получаемого решения и автоматической диагностикой особенностей.
Проведены расчеты, подтверждающие возможность расширения методики для диагностики особенностей при решений систем ОДУ и уравнений в частных производных.
Построенные численные методы применены для интегрирования системы дифференциально-алгебраических уравнений диффузионно-дрейфовой модели полупроводника. Написана программа, позволяющая рассчитывать статические и динамические характеристики полупроводниковых приборов с контролем точности, с различной постановкой граничных условий. В отличие от популярного в настоящий момент алгоритма Шарфеттера-Гуммеля, программа, благодаря разработанным надёжным численным методам позволяет интегрировать систему уравнений модели целиком, без введения упрощающих предположений.
Литература
Sharfetter D.L., Gummel H.K. Large-Signal Analysis of a silicon Read diode oscillator. // IEEE Transactions on electron devices, Vol. 16, NO 1, pp.64-77, Jan. 1969
Shnitnikov A.S., Philatov N.I. Microwave limiter diode performance analyzed by mathematical modeling. // Solid-State Electronics, Vol. 34, NO 1, pp.91-97, 1991.
Самарский А.А. Теория разностных схем. // Наука, 1989, 616 с.
Толстых А.И. Компактные разностные схемы и их применение в задачах аэрогидродинамики. // Наука, 1990, 230 с.
Рогов Б.В., Михайловская М.Н. О сходимости компактных разностных схем // Математическое моделирование, 2008. – Т.20, № 1.- С.99-116.
Shen M.Y., Zhang Z.B., Niu X.L. A new way for constructing high accuracy shock-capturing generalized compact difference schemes // Comput. Methods and Applied Mechanics, 2003, Vol 192, pp.2703-2725.
Radwan Samir F., Comparison of higher-order accurate schames for solving the two-dimensional unsteady Burger’s equation // Journal of Computational and Applied Mathematics, 2005, 174, pp.383-397.
Rosenbrock H.H. Some general implicit processes for the numerical solution of differential equations. // Comput. J. 1963. V.5. №4. P.329330.
Альшина Е.А., Калиткин Н.Н. Вычисление спектров линейных дифференциальных операторов. // ДАН, 2001. – Т.380, № 4. – С.443-447.
Основные результаты опубликованы в работах
Калиткин Н.Н., Корякин П.В. Одномерные и двумерные бикомпактные схемы в слоистых средах // Математическое моделирование, 2009 г., т. 21, № 8, стр. 44-62.
Калиткин Н.Н., Корякин П.В. Бикомпактные схемы и слоистые среды // ДАН, 2007 г., т. 419, №6, с.1-5.
Е.А. Альшина, Н.Н. Калиткин, П.В. Корякин Диагностика особенностей точного решения методом сгущения сеток // ДАН, 2005 г., т. 404, №3, с.1-5.
Е.А. Альшина, Н.Н. Калиткин, П.В. Корякин Диагностика особенностей точного решения при расчетах с контролем точности // ЖВМиМФ, 2005 г., т. 45, №10 с. 1837-1847.
Альшина Е.А., Корякин П.В. Численный метод для режимов с обострениями. // Тезисы докладов II всероссийской конференции памяти А.Ф. Сидорова «Актуальные проблемы прикладной математики и механики», Абрау-Дюрсо, 2004, с. 10-11.
P.V. Koryakin The singularity diagnostics in numerical solving systems of ODE // International congress for mathematicians, Madrid, 2006.
Е.А. Альшина, Н.Н. Калиткин, П.В. Корякин Диагностика особенностей решения при расчетах схемой CROS // Международная конференция студентов, аспирантов и молодых учёных по фундаментальным наукам «Ломоносов-2005», секция физика, сборник тезисов, часть 1, стр. 107-109.
Корякин Павел Владимирович
Моделирование полупроводникового диода
Автореферат диссертации на соискание ученой степени
кандидата физико-математических наук
Отпечатано в Институте математического моделирования РАН
Тираж 100 экз.
|
