Применение методов математического моделирования к исследованию уравнений электрогидродинамики и переноса зарядов в полупроводниках

В настоящее время математическое моделирование получило очень широкое распространение в самых разных областях науки. Огромная сложность наблюдаемых в природе явлений делает невозможным их непосредственное изучение, поэтому вместо интересующего нас явления приходится изучать его модель, в которой отражены все существенные черты данного явления. При описании различных физических явлений часто используются гидродинамические моделиих общность с хорошо изученной моделью газовой динамики позволяет применить разработанные для нее методы к исследованию различных моделей гидродинамики. В частности изучение математических моделей возникающих в релятивистской гидродинамике, магнитной гидродинамике, сверхтекучей жидкости, радиационной гидродинамике, показало, что в ряде случаев систему уравнений описывающую ту или иную модель можно симметризовать (см. 1]), что, в свою очередь, позволяет использовать математически развитую и продвинутую теорию симметрических ¿—гиперболических (по Фридрихсу) систем (хотя на самом деле из-за сложности задач, возникающих на практике, часто требуется дальнейшее нетривиальное развитие этой теории).

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

В последние годы резко возрос интерес к изучению течений сред с объемным зарядом. Дело в том, что заряженные частицы часто появляются в потоках в результате взаимодействия с обтекаемыми телами. Особое значение имеет изучение электризации, возникающей при движении нефти и продуктов ее переработки по трубопроводам (см. 6],[7]). Решение проблемы электризации и снятия зарядов с летательных аппаратов также имеет важное практическое значение. При этом во многих случаях магнитное поле не влияет на движение среды. Изучение движения среды с объемным зарядом в электрическом поле и составляет предмет исследований электрогидродинамики. В обзоре [3] подробно описан широкий круг проблем и приложений электрогидродинамики, а также приведен список работ, в которых формулируются основные предположения, использующиеся при построении математических моделей, описывающих движение сплошной среды с объемным зарядом в электрическом поле. В данной диссертации рассматривается модель, которая описывает движение сплошной среды в электрическом поле. Среда состоит из нейтрального газа и положительно заряженных ионов, причем концентрация ионов много меньше концентрации нейтральных частиц. Пренебрегая вязкостью и теплопроводностью, выпишем уравнения электрогидродинамики. Уравнение неразрывности.

Н + (Иу (ри) = 0, (0.1) где р — плотность сплошной средыи = («1, и2, г/з)* — ее скорость (символ «звездочка» означает транспонирование). Векторное уравнение движения ри)< +сИуП = 0, (0.2) где П — тензор плотности потока импульса с компонентами.

Щ = рщик + р8{к — Р1к (г, к — 1, 2, 3).

Компоненты Р{к максвелловского тензора напряжений Р имеют вид (см. 2]): где р — давление, Е = Еъ, Е^)* — напряженность электрического поля.

Уравнение энергии ре)4 + а^ = (1,Е) — (0.3) где е = во + ||и|2- ео — внутренняя энергия- ^^ = (И7*!, И^, И^)* = /щ (е + руу, v = i/р— удельный объем.

Уравнение сохранения заряда сПуЗ = 0- (0.4) где I — плотность токац — заряд.

Уравнения Максвелла для электрического поля в электрогидродинамическом приближении сНУЕ = 4тгд- (0.5) гс^Е = 0. (0.6) Термодинамические переменные связаны соотношением Гиббса.

Тд, 8 = д. еъ+рдУ (0.7) з — энтропияТ — температура). Из (0.7) вытекают равенства р = -(е0)у = р2{ео)р] Т = (е0)5. (0.8).

Плотность тока Я связана со скоростью и и напряженностью электрического поля Е законом Ома.

Л = д (и + ЬЕ) (0.9) где постоянная Ь > 0 есть так называемая подвижность (см. 12],[3])). С учетом уравнения состояния е0 = еь (/м) — равенств (0.8) и закона (0.9) можно рассматривать систему (0.1)—(0.6) как систему для нахождения компонент векторов и = (р, и*)*, Е и заряда д. В диссертации исследование уравнений электрогидродинамики проведено для случая политропного газа.

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

А (р = -4тгд- (0.10) и в силу этих же соотношений (0.5), (0.6) векторное уравнение (0.2) можно переписать так: ри), + <1М1 = дЕ- (0.11) где П — тензор плотности потока импульса с компонентами Пг£ = рщщ + р5ц-. Система (0.1), (0.11), (0.3) является системой уравнений газовой динамики (с правыми частями), которую можно переписать в недивергентном виде рс2 сИ с2Т йэ, а. .о й д где -— = —- + (и, V) — с2 = (р2(ео)р)р — квадрат скорости звука в газе, а ь с/ с см. 17]). Система (0.12) может быть записана в симметрическом виде.

0.13) 1 где В= и) = с!1а§-(1 /(рс2), 1, р, р, р) — диагональная матрицаВ (к) = В (к)(^ симметрические матрицыЕ = Е (и, Е, д) — вектор правых частей (матрицы и вектор Е без труда могут быть выписаны). Естественно предполагать, что термодинамические величины удовлетворяют неравенствам р > 0, (р2(ео)р)/, > 0. Тогда матрица Вположительно определена, и, следовательно, система (0.13) — симметрическая ¿—гиперболическая (по Фридрихсу)(см. 18],[20]). Итак, с точки зрения теории уравнений с частными производными рассматриваемая модель состоит из симметрической^{гиперболической} системы для р, 5, ад, уравнения для q и уравнения Пуассона для потенциала электрического поля (р.

При исследовании гидродинамических моделей огромное значение имеет исследование сильных разрывов. Такие разрывы встречаются во многих течениях, которые изучаются в механике сплошных сред и представляют большой практический интерес. Сильный разрыв — это поверхность в поле течения среды, при переходе через которую макроскопические параметры среды изменяются скачкообразно. Наиболее часто встречающимся сильным разрывом является ударная волна. Отметим, что сплошная среда после ударной волны подчиняется тому же уравнению состояния, что и среда перед ней. Зона, внутри которой происходит переход от одного физического состояния к другому, имеет конечную, но чрезвычайно малую толщину (порядка нескольких длин свободного пробега). Однако для решения многих задач не требуется подробно знать физические условия внутри переходной зоны. Часто оказывается возможным адекватное описание этого явления путем рассмотрения ударной волны как математического разрыва непрерывности. Итак, математически ударная волна — это поверхность в сплошной среде, которая перемешается с течением времени по частицам среды и на которой параметры среды терпят разрыв первого рода, оставаясь непрерывными с каждой стороны от поверхности. Величины разрывов параметров среды не могут быть произвольными, но с необходимостью удовлетворяют соотношениям, которые получаются из законов сохранения и уравнений Максвелла в интегральной форме.

Если поверхность сильного разрыва задается уравнением где х = (ж^х'), х' = (.г" 2,?з), тогда на этой поверхности должны выполнятся следующие соотношения х)=/(*, х')-я 1 = 0,.

0.14).

Мр] - [рщ] + /х2[рщ] + /Хз[рщ] = о,.

АЫ — [йн] + ЛэРИ + /, 3[п3г-] = о (2 = 1,2,3),.

Ь[ре] - [ТУх] ' д. / [цг3] = 0,.

0.15) (0.16) (0.17).

0.18).

0.19) (0.20) -47ГСГ,.

Л, И = 0 (к — 2,3), где = (Л, К) — ?^ = (Е, 1Ч) — N = —^(-1, Да, Д3)*- 1.

IV/! = /Г+Ж+Ж;

N — нормаль к поверхности (0.14). Кроме этого, мы использовали обычное обозначение.

14 = ^ - ^со! где Гзначение величины ^ справа, при / —> —0, ^— значение величины Е слева, при / —> +0 от поверхности разрыва. Кроме того, для ударной волны ] ф 0, [р] 0, где ] = р{идг — 1) дг), аддг = {и, И), =.

— шл.

При получении соотношений (0.17)-(0.20) предполагалось, что на поверхности (0.14) существует поверхностный заряд, а = стН^х'), и в соответствии с рекомендациями из [3],[12] мы пренебрегли величиной поверхностного тока. Система соотношений (0.15)-(0.20) является замкнутой при заданной величине а.

При численном изучении ударных волн неизбежно возникает вопрос об их структурной устойчивости. Дело в том, что структурно неустойчивая ударная волна при малом возмущении ее фронта разрущается или распадается на несколько ударных волн. Численное исследование таких волн невозможно, поэтому прежде чем переходить к численному исследованию ударных волн необходимо решить вопрос об их структурной устойчивости. Физическая постановка задачи об устойчивости ударной волны состоит в следующем (см. 20]). Пусть сверхзвуковой стационарный поток газа отделен от дозвукового стационарного потока газа поверхностью сильного разрыва — ударной волной (с уравнением х =0). Пусть в начальный момент времени на фронт ударной волны накладывается экспоненциальное возмущение, а набегающий поток остается невозмущенным. Спрашивается, как поведет себя возмущение с течением времени?

Существуют различные подходы к изучению устойчивости сильных разрывов в механике сплошной среды. В газовой динамике первое теоретическое исследование устойчивости плоской ударной волны произвольной интенсивности, распространяющейся в неограниченной среде, было осуществлено в 1954 году С. П. Дьяковым [19]. Суть предложенного С. П. Дьяковым метода состоит в следующем. Сначала система уравнений газовой динамики и условия на сильном разрыве линеаризуются относительно постоянного решения. Затем формулируется некоторая смешанная задача и решение этой линеаризованной задачи ищется в экспоненциальном виде. В результате можно получить так называемое характеристическое уравнение частот. Исследуя это алгебраическое уравнение, Дьяков указал области устойчивых и неустойчивых состояний. В последствии С. В. Иорданский и В. М. Конторо-вич уточнили и дополнили результаты работы Дьякова. В 1967 году Зайдель предложил метод решения смешанной задачи для линеаризованных уравнений, основанный на применении преобразования Лапласа. При этом начальные возмущения могут быть произвольными вдоль одной из пространственных координат (той, которая соответствует нормали к фронту ударной волны). Несколько позже в работе [20] для исследования устойчивости ударных волн был развит метод диссипативных интегралов энергии.

В диссертации использован подход А. М. Блохина (см. 20]), примененный к исследованию задачи об устойчивости ударных волн в газовой динамике. Этот подход заключается в том, что мы формулируем некоторую линейную смешанную задачу, корректность которой изучаем с помощью техники диссипативных интегралов энергии. Если эта смещанная задача окажется корректной, то соответствующий разрыв устойчив и будет существовать как физическая структура. В противном случае разрыв неустойчив и малые возмущения вынуждают его прекратить свое существование даже за незначительное время. Отметим, что в газовой динамике эта методика позволила также перенести результаты исследования линеаризованной задачи на квазилинейные уравнения для изучения гладких течений в окрестности гладкой ударной волны. Для таких течений удалось с помощью построенных диссипативных интегралов энергии доказать теоремы существования и единственности в малом по времени. Также необходимо отметить, что с помощью техники диссипативных интегралов энергии удалось получить ряд важных результатов при исследовании устойчивости ударных волн в релятивистской гидродинамике, гидродинамике сверхтекучей жидкости, в радиационной гидродинамике, в магнитной гидродинамике.

В диссертации рассматривается стационарная электрогидродинамическая ударная волна с малым скачком напряженности электрического поля. Для исследования структурной устойчивости этой волны формулируется линейная смешанная задача, которая получается путем линеаризации уравнений электрогидродинамики (0.4), (0.5), (0.6), (0.12) и соотношений на сильном разрыве относительно рассматриваемого кусочно-постоянного решения. В результате для возмущений неизвестных величин р, 5, и, д получится две симметрические I — гиперболические (по Фридрихсу) системы (слева и справа от разрыва), а нахождение возмущения напряженности электрического поля Е сводится к нахождению возмущения потенциала электрического поля (р. При этом все перечисленные искомые величины связываются граничными условиями. Необходимо также найти Е и О — возмущение ударной волны и поверхностного заряда на ней, для нахождения которых используются два соотношения из граничных условий. Для корректности нашей задачи необходимо проверить, правильно ли задача поставлена по числу граничных условий. Изучение этого вопроса выявило, что возможны несколько вариантов задачи в зависимости от направления тока, возникающего после возмущения ударной волны и поверхностного заряда. Правильно поставленной задача оказывается в тех случаях когда ток течет: 1) через разрыв вниз по потоку без возмущения поверхностного заряда на ударной волне, 2) через разрыв вверх против потока также без возмущения поверхностного заряда на ударной волне, 3) слева и справа в сторону ударной волны, создавая возмущение поверхностного заряда на ней. Все остальные варианты задачи оказывается неправильно поставленными по числу граничных условий, и следовательно, описываемые ими ударные волны являются неэволюционными.

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

Приведем основную идею техники диссипативных интегралов энергии.

Определение 1. Система п уравнений первого порядка с вещественными коэффициентами ои* +? Акихк = 0 — (0.21) к=1 где Аа — квадратные матрицы порядка п, а = 0,3, называется симметрической ¿—гиперболической (по Фридрихсу), если матрицы Аа, а — 0,3, симметрические, а матрица Л о к тому же положительно определена.

Умножим скалярно (0.21) на 2U и получим равенство.

Аои, Щ+Е (А"и, и) Хк=0. к=1.

Предполагая, что |U| = (U, U)1//2 —>¦ 0 при х —> оо и |х'| —> оо и интегрируя это соотношение по области R/j., приходим к уравнению p (t) — J (АМЪ^М- 0, (0.22) Rгде.

I (t)= J (AoU, U) dx. R+.

Для системы (0.21) рассмотрим смешанную задачу с диссипативны-ми граничными условиями.

Определение 2. Граничные условия для симметрической системы (0.21) при х = 0 называются диссипативными, если для всех ненулевых векторов U, удовлетворяющих граничным условиям, выполнено неравенство.

-(AiU, U) >0. (0.23) xi=0.

Если выполнено неравенство (0.23), интеграл (0.22) называется дисси-пативным интегралом энергии. В этом случае (0.22) можно заменить интегральным неравенством т < 7(°) •.

Следовательно, имеет место априорная оценка.

НиМН^^СНиоП^кз);

0.24) где С > 0 — некоторая постоянная, не зависящая от вектора и. Из оценки (0.24) следует, что смешанная задача для симметрических ¿—гиперболических систем корректна поставлена. Однако такие простые рассуждения в нашем случае неприменимы, поскольку граничные условия для основной задачи не являются диссипативными. Чтобы преодолеть это затруднение, в диссертации построена так называемая расширенная система, у которой граничные условия уже диссипативны. Расширенная система — это система для вектора и и производных вектора и, полученная с помощью дифференцирования симметрических систем основной задачи. При построении расширенной системы мы используем то обстоятельство, что из симметрической системы для р, 5, и, д, при х > 0 получается волновое уравнение для р, из которого затем строится симметрическая ¿—гиперболическая система (по Фри-дрихсу), для которой граничные условия уже диссипативны. Задачу, для которой доказана корректность, можно далее исследовать численными методами.

В силу сложности м ат е м ат и ческой модели далеко не всегда можно сделать какие-либо заключения о корректности рассматриваемой задачи. Тем не менее, требуется иметь некоторые представления о поведении решения. В диссертации в рамках предложенной в [22] гидродинамической модели переноса зарядов в полупроводниках рассматривается тестовая задача о п±п-п+ баллистическом диоде (см. [24, 25]), размеры которого порядка 10−5ст, поэтому численный экспериментвозможно единственный на сегодняшний день путь получения информации о решении. Приведем систему уравнений модели переноса зарядов в полупроводниках, взятую из [22].

0.25).

— М + ^-{пу2 + пСТ) + - } = дъ (0.26) о1 ох < тп> д Г, А 2, , 5 М, Я д (2 2 <тпу1 Н—+.

9* 3 тУ.

С2(Т) + ^ + 4 т, 2С2(Т) + = (0.29).

2 т Ь т 2 т> д2: Ф+"(М-п) = 0, (0.30) с где плотность электронов п, скорость электронов г>, напряженность ст, температура электронов Г, тепловой поток ф, электрический потенциал Ф — неизвестные функции от t и х. Величины в правой части равенств (0.26) — (0.29), называемые источниками, имеют вид.

ПУ 67.

92 = — +П — Фх, тр тп.

9 ^ + § С (Т) — с%) д3 = ПУ — Фх — П тп т,.

IV ^ А.

94 =—-Нпу — фх, тта о тп.

1 3, 5 П2(грЛ .5 (ТУ О п /ч ^ ппу ±пуС (Т + -— + — № = ± (§""' + ~пС2(Т) + 2−2- 2 т ш т 2 I т) тч где.

С2(Т).

Кв-Т.

7П т — эффективная масса электрона:

0.26те, И Т0 = 300 К, 0.24те, И Т0 = 77 Кте — масса электрона: те — 9.11 • №~пКд;

Кв ~ постоянная Больцмана:

Кв = 1.38 -10~22 J ¦ К" 1 (1 = КдГ^) 5.

То — температура окружающей среды: Го = 300 К, 77 Кq — заряд электрона: q = 1.6 • 10−19Сс = 11.7 • со, £о ~ диэлектрическая проницаемость: во — 8.85 • 10−12^, Г1 = и гп1 у 1.

Лг (= N0 — ~ плотность легирования;

Оч^а — соответственно плотности донорской и акцепторской примеси. тр = тр (у, Т) — время релаксации импульса, ти) = т-ьи^^Т) — время релаксации энергии, та — та (у:Т) — время релаксации напряженности, тч = тя (у, Т) — время релаксации тепла.

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

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

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

Цель применения стабилизирующего оператора состоит в том, чтобы добиться большего шага по времени, обеспечивающего устойчивый счет. Так как ищется стационарное решение, то точность по времени не важна, а больший шаг по времени уменьшает число итераций до получения стационарного решения, т. е. экономиться время. В результате применения стабилизирующего оператора схема становиться неявной. Это означает, что для нахождения численного решения на новом временном слое требуется решить систему линейных алгебраических уравнений. Размерность данной системы совпадает с числом узлов сетки N. Если матрица линейной системы произвольная, то лучшие методы решения линейных алгебраических уравнений требуют порядка -/V3/2 операций. Но можно расщепить стабилизирующий оператор так, чтобы получение решения схемы на новом временном слое сводилось бы к решению двух линейных алгебраических систем, но с матрицами, позволяющими решать эти системы прогонками, т. е. затрачивать порядка N операций. Неявные схемы требующие порядка N операций для определения численного решения на новом временном слое называются экономичными. Расщепление стабилизирующего оператора выполнено по физическим процессам. Вначале осуществляется только перенос сплошной среды, а затем внутреннее взаимодействие сплошной среды.

Диссертация состоит из введения и трех глав.