Содержание к диссертации
Введение
Глава I. Постановка задачи 26
1.1. Описание элемента термоэлектрической охлаждающей батареи 26
1.2 Система уравнений и граничных условий для расчета теплового процесса в элементе батареи 29
Глава II. Метод решения полной системы уравнений 33
2.1. Описание численного метода 33
2.2. Численное решение в полупроводниковых столбцах N и Р 41
2.3. Аналитическое решение в воздушных и прилегающих к ним областях 47
2.4. Другой способ аналитического решения в воздушной и прилегающим к ней однородных областях 50
2.5. Итерационный алгоритм решения задачи 56
Глава III. Вычислительные ресурсы, используемые в расчетах 59
3.1. Общая структура МВС 1000М 59
3.2. Сетевые решения 50
3.3. Программное обеспечение 61
3.4. Общее программное обеспечение 61
Глава IV. Результаты численного моделирования 61
4.1. Анализ эффективности методики распараллеливания на основе сравнения скорости расчетов 77
4.2. Анализ характеристик производительности термоохлаждающей батареи 78
Заключение
Литература 87
Приложение 93
- Система уравнений и граничных условий для расчета теплового процесса в элементе батареи
- Численное решение в полупроводниковых столбцах N и Р
- Другой способ аналитического решения в воздушной и прилегающим к ней однородных областях
- Анализ характеристик производительности термоохлаждающей батареи
Введение к работе
В данной работе предложен и осуществлен метод решения уравнений теплового баланса, описывающих распределение тепла в термоэлектрическом охладителе с помощью параллельных вычислений.
Термоэлектрическая охлаждающая батарея является многослойной конструкцией, состоящей из материалов с разными коэффициентами теплового расширения, и включает в себя полупроводниковые ветви N- и Р-типа проводимости, припои, диэлектрические теплопереходы. Как правило, батарея состоит из множества электрически последовательно соединенных термопар из полупроводников JV- и Р-типа. При прохождении постоянного тока через термопары происходит поглощение тепла на одних спаях термопар и выделение тепла на других (эффект Пельтъе) [1], а также поглощение тепла внутри термоэлектрических ветвей термопар (эффект Томсона). Остановимся на этих явлениях более подробно.
Классическая теория [47] объясняла явление Пельтье тем, что электроны, переносимые током из одного металла в другой, ускоряются или замедляются под действием внутренней контактной разности потенциалов между металлами. В первом случае кинетическая энергия электронов увеличивается, а затем выделяется в виде тепла. Во втором случае она уменьшается, и убыль энергии пополняется за счет тепловых колебаний атомов второго проводника. В результате происходит охлаждение. С этой точки зрения следовало бы ожидать, что коэффициент Пельтье будет совпадать с внутренней контактной
разностью потенциалов: ^-^,=- -, где Цх и }л2 - энергии Ферми у
контактирующих металлов. На самом деле это неверно. Дело в том, что по классической теории средняя кинетическая энергия теплового движения электронов в обоих контактирующих металлах одинакова. А это неверно
вследствие различного положения уровней Ферми в обоих металлах. Классическое объяснение учитывает только различие потенциальных энергий по разные стороны границы раздела металлов, считая средние кинетические энергии их одинаковыми. Для того чтобы исправить объяснение, надо изменение потенциальной энергии при переносе электрона из одного металла в другой заменить изменением полной энергии. Исправленное таким образом объяснение, разумеется, справедливо не только для металлов, но и для полупроводников с электронной проводимостью.
Совершенно аналогичное объяснение можно привести и для того случая, когда в контакте находятся два полупроводника с дырочной проводимостью. Через границу раздела переходят, конечно, электроны. По одну сторону границы происходит рождение пар электрон-дырка, по другую - рекомбинация электронов с дырками. Один из этих процессов сопровождается выделением, другой - поглощением энергии. От соотношения между выделяющейся и поглощающейся энергией зависит знак коэффициента Пельтье.
Эффект Пельтье, как и все термоэлектрические явления, выражен особенно сильно в цепях, составленных из электронных и дырочных полупроводников. Рассмотрим контакт таких полупроводников. Допустим, что электрическое поле имеет такое направление, что ток идет от дырочного полупроводника к электронному. Тогда электроны в электронном полупроводнике и дырки в дырочном будут двигаться навстречу друг другу. Электрон из свободной зоны электронного проводника после прохождения через границу раздела попадает в заполненную зону дырочного полупроводника и там занимает место дырки. В результате такой рекомбинации освобождается энергия, которая и выделяется в контакте в виде тепла. Рассмотри теперь случай, когда ток проходит через границу раздела от электронного полупроводника к дырочному, Тогда
электроны в электронном и дырки в дырочном полупроводниках будут двигаться в противоположные стороны. Дырки, уходящие от границы раздела, будут пополняться в результате образования новых пар при переходах электронов из заполненной зоны дырочного полупроводника в свободную. На образование таких пар требуется энергия, которая поставляется тепловыми колебаниями атомов решетки. Электроны и дырки, образующиеся при рождении таких пар, увлекаются в противоположные стороны электрическим полем. Поэтому пока через контакт идет ток, непрерывно происходит рождение новых пар. В результате в контакте тепло будет поглощаться. Таким образом, если ток идет от дырочного полупроводника к электронному, то тепло Пельтье выделяется. При обратном направлении тока оно поглощается. Это будет наглядно продемонстрировано в наших дальнейших расчетах.
А.Ф. Иоффе предложил [1] использовать явление Пельтье в полупроводниках для создания охлаждающих устройств. Как будет показано из анализа эффективности работы ТЭОБ, термоэлектрический метод охлаждения обладает рядом преимуществ по сравнению с другими методами охлаждения; главное из этих преимуществ — малость геометрических размеров термоэлектрического охладителя.
Дополнительным термоэлектрическим эффектом в полупроводнике при прохождении через него электрического тока является поглощение или выделение так называемого «тепла Томсона». Суть данного явления в том, что различные участки термопары нагреты неодинаково, а потому их физические состояния также неодинаковы. Неравномерно нагретый полупроводник должен вести себя как система находящихся в контакте физически разнородных участков. На этом основании Томсон пришел к заключению и подтвердил его экспериментально, что не границах таких участков должно происходить выделение или поглощение тепла Пельтье, Такое тепло получило название тепла Томсона, а само явление — явления Томсона.
-^ E
gr-аіі Т
*" J(oxna»№KH) I
и - иолрироводнык
а)
р - иол^яроводиик'
б)
Рис.1
С точки зрения электронной теории явления Томсона объясняется просто [47]. Рассмотрим полу проводник с электронной проводимостью. Пусть Тх >Т2, т.е. градиент температуры направлен отточки 2 к точке 1 (рис. 1а). Из-за диффузии концентрация электронов в точке 1 сделается меньше, чем в точке 2. Возникнет электрическое поле Е, направленное от 1 к 2, т.е. против градиента температуры. Если по проводнику течет ток в направлении grad Т (т.е. электроны двигаются в направлении поля Е), то поле Е будет замедлять электроны, а участок полупроводника 1-2 станет охлаждаться. Если же ток течет в обратном направлении, то произойдет нагревание участка 1-2. В дырочном полупроводнике соотношения будут обратными (рис. 16). Явление выглядит так, как если бы на обычные поток тепла, вызванный теплопроводностью, накладывался дополнительный поток тепла, связанный с прохождением электрического тока. В дырочных полупроводниках дополнительный поток тепла направлен в ту же сторону, куда течет электрический ток; в электронных направления электрического тока и тепла противоположны. Эффект Томсона считается положительным, если электрический ток, текущий в направлении градиента температуры, вызывает
нагревание полупроводника, и отрицательным, если при том же направлении он охлаждает полупроводник.
Термоэлектрическая охлаждающая батарея, таким образом, представляет собой последовательность термопар, соединенных между собой таким образом, чтобы эффекты Пельтье на спаях и Томсона внутри полупроводниковых ветвей были отрицательными и усиливали друг друга, с целью охлаждения исследуемого объекта.
Термоэлектрические охладители широко применяются в различных областях техники, в том числе в таких устройствах, где тенденция к миниатюризации налагает существенные ограничения как на размеры самих охладителей, так и на возможности отвода тепла от «горячих» спаев термоэлементов (см.[2]).
При существенно ограниченном теплоотводе задачу расчета и проектирования термоэлектрических охладителей необходимо решать в комплексе с задачей о теплообмене с окружающей средой. Кроме того, значительное влияние на тепловое поле в микрохолодильнике может оказать тепловыделение на коммутациях, припоях, подводящих проводниках и т.д. В связи с этим аналитическое рассмотрение задачи оказывается неприемлемым ввиду чрезвычайной громоздкости получаемых при расчете выражений.
Постараемся кратко проиллюстрировать основной существующий подход к моделированию полупроводниковых термоохладителей.
Общие традиционные принципы моделирования полупроводниковых приборов приведены в работе [23]. В работе [23] дается обзор достижений в
области моделирования полупроводниковых приборов. Вначале приводятся основные уравнения поля с соответствующими граничными условиями. Затем описываются эмпирические модели физических механизмов, определяющих работу приборов, например, подвижности, лавинного умножения носителей, сужения запрещенной зоны. Далее кратко излагаются различные численные модели, большинство из которых разработано за последние два десятилетия, и обсуждаются методы дискретизации и решения уравнений моделей. Приводится ряд замечаний по поводу взаимосвязи между методами конечных разностей и конечных элементов. Рассматриваются упрощенные численные модели, и подчеркивается их полезность для определенных приложений. На ряде характерных примеров наглядно показываются возможности численного моделирования полупроводниковых приборов.
Наиболее распространенным подходом к описанию электрофизических процессов в полупроводниковых элементах является подход, основанный на использовании одночастичных функций распределения носителей заряда / - го сорта f,{p,7j) по импульсам р, координатам ? и времени /, где индекс / определяет принадлежность к той или иной долине или подзоне зоны проводимости и валентной зоны,
В начале работы [23] приводятся основные уравнения физических процессов в полупроводниковых приборах и граничные условия. Можно показать, что для полупроводниковых приборов с линейными размерами в несколько микрометров и меньше и в диапазоне частот до 1019 герц справедлива квазистатическая аппроксимация уравнений электрического поля. К ним относятся уравнения непрерывности для плотности электрического тока J:
v,~22 (1)
и уравнение Пуассона для электростатического потенциала у/:
Д^ = - (2)
где - диэлектрическая проницаемость.
Плотность пространственного заряда р складывается из плотностей зарядов
подвижных носителей двух типов (электронов п и дырок р) и
результирующей примесной концентрации Л":
p = q{p-n+N) (3)
где д - (положительный) электронный заряд, а результирующая примесная
концентрация равна сумме полных концентраций ионизированных доноров ND
и акцепторов NA:
N = Nn~NA. (4)
Для полупроводников полный перенос зарядов складывается из переноса
электронов в зоне проводимости и переноса дырок в валентной зоне:
J = J„+JP (5)
Такое разделение требует введения отдельных уравнений непрерывности для каждого типа носителей
-lv.y.=-*; (6)
ot д
3%1 dt q
^ + -^-/, =~*. (7)
где R - результирующая скорость рекомбинации. В изотермических условиях
плотности электронного и дырочного токов пропорциональны градиентам
соответствующих электрохимических потенциалов, называемых
квазиуровнями Ферми:
Зп=-<№п*ЯФп', (8)
J,=-qMppV4p, (9)
где jl - подвижность.
На идеальном контакте граничными условиями для электростатического потенциала и для концентрации носителей обоих типов являются условия Дирихле:
=о+Уар,*\ (10)
п = п0\ (11)
Р = Ро> (12)
где Vappi - приложенное к контакту напряжение; ц/й, п0 и р0 - значения
соответствующих переменных при нейтральности пространственного заряда и в равновесных условиях,
В случае неидеальных контактов, например, контактов Шоттки, граничные условия имеют вид (см. [23]):
n-Jn = ~qv„(n -щ); (13)
n-Jp=
где vn, v - скорости термоэлектрической рекомбинации у контакта Шоттки;
Я - единичный вектор, перпендикулярный границе раздела «металл — полупроводник». При бесконечно больших скоростях рекомбинации (идеальный контакт) эти уравнения переходят в (11) и (12).
Формулы (1)-(14) являются основополагающими выражениями, описывающими механизм переноса электронов и дырок в полупроводниковых приборах. Для решения систем таких уравнений, дополненных другими, частными, выражениями, используются методы дискретизации, сводящие исходную систему дифференциальных уравнений к алгебраической задаче в пространстве с конечной, хотя и достаточно большой размерностью. Из методов дискретизации, то есть вывода разностных уравнений, широкое распространение получили метод конечных элементов и метод конечных разностей. Метод конечных разностей основан на локальной аппроксимации
дифференциального оператора некоторым разностным оператором, а в методе конечных элементов искомое решение глобально аппроксимируется набором функций формы, которые служат пробными функциями. При этом в процессе аппроксимации применяется разностный метод (Галеркина) или вариационный метод (Ритца) [23]. Обоим методам присуще разбиение всего объема прибора на подобласти и использование пробных функций.
Далее в работе [23] приводится иллюстрация применения метода конечных разностей и конечных элементов для дискретизации уравнения Пуассона, которое в интегральном представлении имеет вид:
fr-dS = --\pdv (15)
gv Є У
для любой подобласти v.
Отмечается также важная роль, которую играет согласование приближенных значений с физическим смыслом задачи. Это особенно важно в случаях, когда свойства решений определяются нелинейностями. Иными словами, применяемые пробные функции должны, прежде всего, отражать соответствующие физические процессы, независимо от того, служат ли они одновременно функциями формы или нет. Поэтому целесообразно выбирать в качестве пробных функций приближенные аналитические значения исследуемой задачи.
В качестве иллюстрации в работе [23] приводится анализ уравнения непрерывности для дырок (7), который приводит затем к аппроксимации Шарфеттера-Гуммеля [20] для плотности дырочного тока:
Фр,і*і Фр,1
';+ї ( \ W,^-V, ет -ет
е ' -е
Ч 2 -V] -*(. — —
Также очень удобна аппроксимация:
1Ґ Л
= --<*/+]-*mM*J
Данная разностная схема обычно имеет значительно более широкие пределы применимости, чем другие, более традиционные схемы, в которых плотности тока определяются из полиномиальных пробных функций для квазиуровня Ферми или для концентрации носителей. Если принять \ц/ы -у,|« VT и
|^,*~^j+i|<
Разностную схему Гуммеля можно обобщить и на многомерный случай. В этой ситуации составляющие плотности тока вдоль линий сетки в промежутках между узлами принимаются постоянными и аппроксимируются выражениями вида (16). Именно на базе схемы Гуммеля были выполнены первые работы по точному двумерному численному моделированию биполярных приборов. Такая схема пригодна как для прямоугольных, так и для треугольных сеток.
В литературе описаны и другие примеры физического подхода к дискретизации основных уравнений; например, в работе [42] на базе физических представлений получена точная аппроксимация механизма рекомбинации Шокли-Рида-Холла, а в работе [44] первый интеграл из одномерного уравнения Пуассона включен в разностную схему для приповерхностных областей с окисным покрытием.
В работе [22] рассматривается распространение разностной схемы Шарфеттера - Гуммеля на уравнение энергетического баланса в одномерном приближении. В работе выводятся сеточные значения электронного тока JM и
потока тепла S
Jm=~
'
хм xt
xi+l xt
МТ)1+1 схр[(ам +1)1п(7;+1/Г,)]-Мг)1
laM+lMTM/Tt)]-l
где aM = - g(y ,+1 - y,)1 кB(Tl+l-Tt); tcB - постоянная Больцмана;
Ч
хм xt
&]{уі+і -,) , Ькв(Тм-Т)
хм хі
„ [пМЩ1] ехр[(«м + Ад-1)\п(Т^Т,)]-[пМ{т)т1
(ам + ^8'^n{TiJTi)\-l
И И
где 6 = <>(Т)=(тє7)/{тє},
Д = Д(Г) =
{«)
(e).
«энергозависимое» время релаксации отдельного столкновения.
Работа [23] хорошо дополняется работой [18], которая также производит общий обзор математического моделирования как метода исследования электронных процессов в элементах интегральных схем (ИС), определяющих электрические характеристики и их особенности. Основное внимание в работе [18] уделяется математическому моделированию элементов ИС с субмикронными размерами технологических неоднородностей. Именно такие элементы, по мнению авторов статьи, являются основой перспективных полупроводниковых приборов, обеспечивая последним высокую степень интеграции и быстродействие.
Отметим, что на субмикронных элементах наблюдается сильная неравновесность их электронно-дырочной плазмы, обусловленной
соизмеримостью размеров технологических неоднородностей / и длин, характеризующих релаксационные процессы в ней. К таким характерным длинам относятся, например, длина энергетической релаксации носителей заряда Ле, длина релаксации импульса Я, длина парных столкновений Л. Наиболее сильно эффекты неравновесности выражены в элементах на основе материалов с хорошими динамическими параметрами носителей заряда, такими, как GaAs, InP и др. Однако, и в субмикронных элементах на основе кремния могут существовать области, в которых эффекты неравновесности весьма значительны. Важно также, что неравновсность электронно-дырочной плазмы в субмикронных элементах ИС обусловливает и нелокальность связи между напряженностью электрического поля и его параметрами.
Настоящая диссертационная работа не учитывает, однако, субмикронных неоднородностей полупроводниковых элементов, поскольку, во-первых, рассчитываемая (и реально применяющаяся в технике) ТЭОБ имеет размеры, намного превосходящие субмикронные, а во-вторых, добавление дополнительных «уточняющих» уравнений в исходную систему, хотя и повышает точность модели, но вместе с тем значительно увеличивает машинное время компьютерных расчетов, которое и без «субмикронных поправок» занимает в каждом случае несколько часов работы ЭВМ.
В работе [17] предложен метод численного решения системы уравнений физики полупроводников, описывающей поведение неравновесной электронной плазмы в полупроводниковых структурах. Работа подчеркивает, что ставшая уже классической в численном моделировании элементов интегральных схем («ИС») дрейфово-диффузионная модель [34], область применимости которой ограничена элементами с размерами, превышающими
длины релаксации импульса и энергии носителей заряда, не позволяет
учитывать неравновесность и нелокальность электронно-дырочной плазмы
(«ЭДП»). Несомненно, что для достаточно широкого класса приборов
микроэлектроники эта модель дает приемлемые результаты, особенно если нас
интересуют только выходные характеристики элементов. Однако,
«внутренние» физические процессы, предсказываемые диффузионно-дрейфовой моделью и температурными моделями, принципиально различны. По-видимому, этим можно объяснить тот факт, что с начала 1980-х годов стали интенсивно развиваться подходы, связанные с использованием температурных моделей. Наиболее общее разделение таких моделей можно провести на гидродинамические и квазигидродинамические, которые справедливы в различных, вообще говоря, непересекающихся физических ситуациях. Квазигидродинамическая модель представляет собой сильно нелинейную систему параболических уравнений, решаемых совместно с уравнением Пуассона для потенциала электрического поля носителей заряда.
Данная работа [17] посвящена разработке методов численного решения квазигидродинамической модели применительно к полевым транзисторам с затвором Шоттки (ПТШ) на арсениде галлия. Рассчитывались кремниевые и арсенидгаллиевые ПТШ с субмикронными размерами. Вычислялось значение тока на стоке ПТШ, скорости дрейфа и эффективная температура вдоль канала ПТШ. Показано, что в значительной области скорость дрейфа электронов превышает и скорость насыщения, и максимальную равновесную скорость дрейфа; в области между затвором и стоком эффективная температура электронов существенно отличается от температуры решетки, что связано, очевидно, с наличием в этой области сильного электрического поля, Как показал анализ, распределение эффективной температуры сильно отличается от распределения, полученного из простого уравнения баланса мощности. Это свидетельствует о существенности процессов электронной теплопроводности и
конвективного переноса электронной энергии, что обусловливает ограниченность применения при моделировании субмикронных арсенидгаллиевых и кремниевых ПТШ дрейфово-диффузионного приближения.
В работах [2-4] были представлены попытки расчета распределения температуры в различных элементах термоохладительных батарей в каждой точке расчетной сетки. Были получены численные решения системы уравнений теплового баланса. Вместе с тем, в работе [2], в методе численного решения данной системы не принималась в расчет двумерность модельной задачи, т.е. стационарная температура в различных слоях полупроводниковых каскадных термоохладителей полагалась одинаковой вдоль каждого горизонтального слоя. Кроме того, в качестве необходимых исходных данных для решения задачи, предложенного в работе [2], фигурируют величины, которые зачастую нелегко установить с требуемой точностью, например: эффективный коэффициент теплообмена охлаждаемого объекта, интегральный коэффициент теплообмена теплоотвода со средой. В модельной задаче налагаются ограничения на форму каскадного термоохладителя, что сужает область его применимости.
При решении главной задачи - определения температурных полей внутри термоохладителя - возникает вопрос о точности математической модели, т.е. об исходной системе уравнений, описывающих тепловые процессы в термоохладителе, и о задании начальных данных. Объективно говоря, в качестве начальных данных желательно использовать разность потенциалов, приложенную к контактам термоохладителя, и затем, используя одну из вышеуказанных моделей, переходить к вычислению протекающего электрического тока и к сопутствующим тепловым эффектам. Однако, чрезвычайная сложность такой модели в ее представлении как компьютерной расчетной программы затрудняет использование эффективных методов
стыковки и распараллеливания задачи на многопроцессорной ЭВМ, чрезвычайно увеличивает машинное время расчетов и сужает допустимые размеры расчетных областей.
В связи с вышеизложенным было решено в качестве исходных данных рассматривать силу электрического тока, протекающего через термоохладитель, а основным уравнением, описывающим тепловые процессы в нем, считать стационарное уравнение теплового баланса, которое в дифференциальной форме имеет вид в нашем случае: div(*Vr)+ pf - TjVa + QPeUier = О
(пояснения см, в разделе 2 Главы I).
Отметим, что в настоящем уравнении учтены поток тепла в результате теплообмена, джоулево тепловыделение, а также тепловые эффекты Пельтье и Томсона, возникающие при прохождении тока через спаи различных материалов и через полупроводник, в котором присутствует отличный от нуля градиент температуры.
В работе [11] был произведен анализ применимости метода приближенной
факторизации (по сути тот же, что и в настоящей диссертационной работе) для
решения уравнения Пуассона, описывающего распределение
электростатического потенциала в двумерных областях типа «металл — диэлектрик - полупроводник» со сложной геометрией и глубокими профилями металла, окисла и примеси.
Метод приближенной факторизации широко используется при расчетах потенциальных течений сжимаемого газа [27]-[29]. При этом для задач с сильной нелинейностью в правой части факторизованный оператор можно выбрать таким образом, что вторые производные разделены, и при реализации алгоритма каждая его часть обращается изолированно на своем шаге итерационного процесса. Получаемый таким образом факторизованный метод
известен в литературе как метод ПФ1 (см. [27]). При расчете потенциальных течений сходимость такой схемы зависит от типа оператора и сильно ухудшается в области сверхзвуковых скоростей, когда тип уравнения становится гиперболическим. Необходимость в специальных приемах для обеспечения устойчивости процесса ограничивает область применения ПФ1 для решения такого класса задач. Однако, при расчете электрических потенциальных полей, когда уравнение имеет строго эллиптический тип, использование данного метода становится перспективным (см. [11]).
Одной из основных работ, оказавших влияние на выбор методики решения рассматриваемой модельной задачи, является работа [10]. В данной работе предлагается методика построения оптимальных итерационных схем приближенной факторизации для решения двумерных уравнений, описывающих установившиеся трансзвуковые потенциальные течения. Под оптимальными схемами понимаются итерационные факторизованные схемы с наиболее высокой скоростью сходимости. Показано, что при построении оптимальных факторизованных схем для уравнения потенциала необходимо учитывать преобразование координат и правильно распределять производные функций преобразования между двумя сомножителями факторизованного оператора. В работе находятся оптимальные значения итерационного параметра, определяющего устойчивость схемы, и описывается методика выбора оптимальной последовательности значений итерационного параметра, определяющего скорость сходимости итерационной схемы. Рассматривается задача решения разностных уравнений следующего типа:
W) = (ASXSX + BtiZSz)fi = Q, (17)
где 5Х и SZ разности для интерполирования вперед (например,
ЗХФі - Фм ~Фі), ь ёХ и $2 - разности для интерполирования назад (например,
дХфі =Фі~ фіЛ). Коэффициенты Л и В в общем случае предполагаются
переменными.
Итерационную схему решения разностных уравнений (17) можно записать в
виде:
где А = ф<п*і] ~фіп) - «поправка», вычисляемая на n+7-й итерации, а с -итерационный параметр.
В схемах приближенной факторизации оператор N представляется в виде произведения двух сомножителей: например, простейшая форма записи схемы приближенной факторизации 1-го типа [27] имеет вид:
(а - А ЗХ<5х\а - В <5Z^z)i = ааі(фм ),
где а - итерационный параметр, от выбора значений которого зависит
скорость сходимости.
На практике решение приведенного выше уравнения сводится к последовательному решению уравнений: (а-ЛвХ$х)р = асгі{ф(")), (a-BZSZ.)\ = F.
Такой итерационный процесс быстро сходится при оптимальной последовательности значений а. При расчете течения с помощью схемы переменных направлений возмущения быстро распространяются по всему полю течения, и поэтому ее скорость сходимости может быть очень большой. Однако, для реализации этой возможности нужно правильно разложить оператор N на сомножители и выбрать оптимальный набор значений итерационного параметра а.
В работе [7] были теоретически найдены условия сходимости итерационного метода, позволяющего сводить решение краевой задачи для уравнения
Пуассона в области, составленной из нескольких цилиндрических областей, к решению на каждом итеративном шаге краевых задач отдельно в каждой из составляющих областей. Как будет показано ниже, ряд приемов данной работы [7] может быть применен к расчету тепловых процессов в рассматриваемой нами модельной задаче.
В работе [8] исследуется и обосновывается допускающий полное распараллеливание итерационный метод с разделением на подобласти (без перекрытия) решения краевых задач для уравнения Аи+ц2и = /, где А -положительный эллиптический оператор второго порядка, р - большой параметр. Рассматривается случай цилиндрической области с липшицевым сечением, на боковой поверхности задаются смешанные условия. Разбиение производится по гиперплоскостям, параллельным торцам. Метод сходится со
скоростью геометрической прогрессии со знаменателем порядка е м , где / — минимальная ширина слоев разбиения. В обсуждаемой здесь работе [8] рассматривается следующий случай: пусть область Q. -цилиндр Q = Ax(a,b) в и-мерном евклидовом пространстве 91" точек
x={xt х„), где Л- ограниченная область с липшицевой границей в
пространстве їИя"1точек х =(х„ ,*„_,), а<Ь. Граница Г такой области Q
также является липшицевой (см. [31] и библиографию к ней). Исследуемый метод состоит в разрезании области fi на т>\ подобластей Qy,
j=I,2,....,m, представляющих собой слои
(*) їїj = Лх(Я>1,at),j=J,2 т,
где а = а0 < <ат =Ь. Обозначим через 1} =(а, -я;_,) ширины слоев п,, через
Л(=Лх{й(}, і=1,2, ,тЛ,
Система уравнений и граничных условий для расчета теплового процесса в элементе батареи
В качестве потенциального добавления к настоящей работе, исследующей специальный численный метод для расчета тепловых эффектов, можно привести результаты, описываемые в работе [19]. Автор работы справедливо отмечает, что электрический ток в полупроводниках при однородной температуре решетки TJj, но неоднородным распределением электрического поля (например, при наличии барьерного слоя) состоит из слагаемых, отвечающих за проводимость, диффузию и термическую диффузию. Первые двое слагаемых содержат в себе подвижность и коэффициент диффузии электронов, которые являются функциями электронной температуры Г, или, в общем случае, зависят от определенных усредненных величин на неравновесной, зависящей от величины электрического поля, функции распределения энергии электронов. Третье слагаемое соответствует градиенту электронных температур и аналогично общеизвестной термической диффузии газа в условиях градиента температур, В общепринятой теории, которая пренебрегает нагревом или охлаждением электронов, подвижность и коэффициент диффузии являются материальными константами, а термическая диффузия отсутствует, В противоположность случаю равномерных электрических полей, Т не является функцией единственно локального поля; она также зависит от электрического тока и может быть определена лишь при одновременном решении уравнений для потока носителей зарядов и сохранения энергии с граничными условиями на конкретной структуре. В качестве примера, работа [11] рассматривает выпрямитель, основанный на контакте металла и полупроводника с носителями заряда одного типа; детальный анализ включает в себя также обсуждение эффекта Пельтье. В области барьера Т больше, чем То (т.е. эффект «горячих электронов») для обратного смещения, но меньше, чем То, (т.е. «холодные электроны») для прямого смещения. В работе [19], в частности, исследуется с помощью компьютерного моделирования поведение электронов при прохождении барьера Шоттки, учитывая рассеяние электронов на акустических фононах. Расчет всей термоохладительной батареи требует большого объема вычислительных ресурсов и памяти ЭВМ ввиду ее значительной длины. Поскольку описываемая батарея состоит из большого количества одинаковых элементов, имеет смысл исследовать возможность одновременного расчета этих элементов на многопроцессорной ЭВМ.
В связи с этим возникает задача о разделении всей рассматриваемой области на подобласти (без перекрытия). Необходимо выбрать границы подобластей и метод стыковки решения.
Метод, изложенный в данной работе, был апробирован автором на модельном элементе аналогичного типа на однопроцессорной ПЭВМ [61]. Настоящая же работа описывает распространение метода на методику параллельных, многопроцессорных расчетов.
Следует отметить, что данной тематике при рассмотрении краевых задач для эллиптических уравнений посвящено большое количество работ. Автор хотел бы особо отметить работы [6-8], под влиянием которых выполнялось данное исследование. В настоящей работе предложен итерационный метод разделения подобластей со стыковкой на границах решения, учитывающий специфику и особенности задачи. Целью работы является разработка вычислительной методики, позволяющей моделировать тепловые процессы в термоэлектрической батарее с использованием многопроцессорного компьютера с распараллеливанием вычислений в различных отдельных элементах.
Предлагаемый метод был успешно реализован при моделировании данной задачи на многопроцессорном вычислительном комплексе «Суперкомпьютер МВС 1000М». Результаты расчетов и возникшие при этом особенности также приведены в настоящей работе.
Работа состоит из Введения, четырех глав, Заключения, Списка литературы и Приложения. Первая глава посвящена постановке задачи распределения тепла в термоэлектрической охлаждающей батарее - полупроводниковой структуре, в двумерном случае. Приведены модели используемых электрофизических параметров, рассчитано распределение электрического тока в различных частях термоэлектрической батареи, учтены эффекты Пельтье и Томсона. Рассмотрены различные типы граничных условий.
Численное решение в полупроводниковых столбцах N и Р
Она состоит из шести типов различных частей, разделенных вертикальными прямыми 1,2,3,4,5, ...,п. Области (1-2), (2-3), (3-4), (4-5),...,(п-1,п) многократно повторяются в реальной термобатарее.
Области (1-2), (3-4),...,(п-1,п) состоят из полупроводниковых столбцов различной проводимости и участков припоя, контактов и алюминиевой подложки с тонким слоем диэлектрика. Поскольку в полупроводниковых столбцах уравнение (1.4) нелинейно, то решение в этих областях можно получить только численно. Назовем эти области в дальнейшем полупроводниковыми.
Области (2-3), (4-5),...,(п-2,п-1), изображенные на рис.6 как в-в (воздух -внутренний), состоят из воздушного зазора, а также участков припоя, контакта и алюминия с диэлектриком. Такие области назовем областями аналитического решения.
Области, изображенные на рис.6 как в-л (воздух — левый) и в-п (воздух -правый), являются левой и правой соответственно воздушными краевыми областями и искусственно вводятся для задания температуры на левой и правой границе этих областей. Предполагается, что эти области достаточно протяженны, чтобы процессы теплообмена, проходящие в батарее, не влияли на температуру окружающей среды.
Таким образом, используется итерационный метод решения смешанных краевых задач с разделением на подобласти (без перекрытия).
Стыковка решений проводится по вертикальным участкам границ подобластей 1,2,3,4,5,...,11 и подобных им, для которых приводятся специальные граничные условия и алгоритм стыковки.
Поскольку расчет задачи был осуществлен на многопроцессорном компьютере, то схожая последовательность вычислительных действий производилась параллельно (одновременно) несколькими процессорами. Каждому процессору соответствовал определенный участок расчетной области (назовем его в дальнейшем здесь расчетным блоком); схема разделения показана на рис.6 (области, выделенные снизу дугами, соответствуют различным расчетным блокам).
В первоначальный момент во всей расчетной области задается начальное распределение температуры TQ(x,y) (более подробно о задании начального распределения расскажем в п. 1а).
В полупроводниковых столбцах решается уравнение (1.4) со смешанными, выбранными специальным образом граничными условиями. Как уже говорилось, расчет в электронных и дырочных полупроводниковых столбцах (области N и Р на рис.6) проводился одновременно, параллельно для всех частей.
Решение для температурных полей в подобластях в-л и в-п, а также во всей промежуточной подобласти в-в ищется аналитически, так как в этих областях уравнения теплового баланса являются линейными. В воздушных областях, в алюминиевой подложке с диэлектриком правая часть уравнения теплового баланса равна нулю, а коэффициенты теплопроводности - постоянные величины. В области припоя и контакта правая часть определяется заданной величиной плотности тока, и уравнение теплового баланса (1.4) в этих областях также линейно. При этом на вертикальных границах 2,3,4,5, ...,п-1. промежуточных областей (в-в) и на границах 1,п для воздушных областей в-л и в-п задаются значения температуры, полученные из численного решения уравнения (1.1) в областях (1-2), (3-4),,,,,(п-1,п). Поскольку производные температуры по вертикали на этих границах имеют разрыв из-за того, что в областях (1-2),(3-4),...,(п-1,п) коэффициенты теплопроводности рвутся, промежуточные области в-в, области в-л и в-п делятся на семь подобластей, границы которых отмечены штриховыми горизонтальными прямыми. В каждой из этих подобластей строится свое аналитическое решение для Т. При этом в распределении Г выделяется линейная часть, а нелинейный остаток раскладывается по гармоникам ряда Фурье.
В первоначальных вариантах расчетных схем мы не использовали аналитический поиск решения в воздушных подобластях «в-л», «в-п» и «в-в». Единый метод (описанный ниже в п.2) использовался для поиска решения во всех без исключения подобластях модельной задачи (происходил «сквозной счет» по полупроводниковым, проводниковым, и воздушным областям).. Однако, при этом мы столкнулись с трудностями, заключавшимися, прежде всего, в крайне медленной сходимости расчетного метода. Низкая скорость сходимости была обусловлена существенной разницей коэффициентов теплопроводности (до 4-х порядков) между токопроводящими областями и воздухом. Было очевидно, что при такой медленной «подгонке» тепловых потоков расчетный метод не может считаться эффективным. В связи с этим и было предложено дополнить расчетный метод аналитическими разложениями в воздушных (вакуумных) и прилегающих к ним сверху и снизу подобластях.
Поскольку диэлектрический слой между алюминиевой подложкой и припоем весьма тонок (в 10 раз тоньше каждого из них), то мы рассматриваем алюминиевую подложку вместе с диэлектриком как однородный слой с эффективной теплопроводностью [6]
Другой способ аналитического решения в воздушной и прилегающим к ней однородных областях
В работе использовались экспериментальные значения ряда физических величин [13]: теплопроводность воздуха лгл = 2.41x10 Дж/см -с-К; теплопроводность алюминия к =2.302 Дж/см -с-К; теплопроводность диэлектрика л"д =1.67хЮ 3Дж/см -с-К; теплопроводность припоя кп = 4.51x10" Дж/см -с-К; теплопроводность шины-контакта к,0„ =3,895 Дж/см -с-К; удельное сопротивление припоя ри =11.5 :10"6Омсм; удельное сопротивление контакта р% = 1.56x10"" Ом-см; площадь поперечного сечения припоя 5л=і.4хіо3 см2; площадь поперечного сечения контакта 5,=4.2хЮ3 см ; площадь поперечного сечения термоэлектрических ветвей Snp =I.96xI0 2 см2. Кроме того, была проведена дополнительная серия расчетов, при которой рассматривался «случай вакуума», то есть теплопроводность воздуха была условным образом уменьшена до кь = 2.41x10" Дж/см -с-К.
Эффективная теплопроводность алюминиево-диэлектрического слоя в зависимости от толщины диэлектрика равна значениям, приведенным ниже: Все эти данные, а также теплопроводность, удельное сопротивление и термо-ЭДС термоэлектрических ветвей приведены в [13]. Из рассмотрения рис.16 можно оценить величину {daf&r\p «8x10" Дж/Кл-К. При проведенных расчетах определяется влияние силы электрического тока на источник тепла и термо-ЭДС в составляющих элемента, а также на тепло Пельтье на спаях полупроводниковых столбцов. Изучается зависимость распределения температуры от задания внешнего теплового потока Q. Расчеты проводились при следующих допущениях: температура окружающей среды вблизи "холодной" поверхности термобатареи равнялась 260 К, температура вблизи "горячей" поверхности 290 К. Вдоль левой границы области (в-п) „ =290 К, вдоль правой границы области (в-п) - аналогичным образом: п =290 К. Расчеты проводились для значений тока, равных 7=0.1 А, 0.3А, 1А, ЗА, 4А, 5А. При этом выяснилось, что задача остается физически оправданной лишь при известных значениях теплового потока на «холодной» границе термобатареи, причем этот изначально задаваемый поток различен для разных значений тока J (см. также рис.4 (профиль задаваемого потока) и комментарий на стр. 32).
Расчетная программа, составленная в ходе настоящей работы, позволяет произвольно изменять чрезвычайно большое количество физических параметров модельной задачи. Например, можно изменять размеры всех составляющих элементов термоэлектрической батареи, их физические свойства (т.е. изменять материал, из которого изготовлены элементы), силу тока, протекающего через ТЭОБ, граничные условия (поток тепла на верхней поверхности и температуру на нижней) и пр. В связи с этим после проведения большого количества различных расчетов можно, по всей видимости, подобрать ряд оптимальных параметров. Однако, основная цель настоящей работы заключается не столько в конструировании термоэлектрической охлаждающей батареи, сколько в разработке численного метода, который позволит рассчитывать процессы в конструкциях, состоящих из большого количества одинаковых (или, во всяком случае, схожих между собой) элементов. Данный приоритет целей заставил нас ограничить число изменяемых параметров модельной задачи — до пяти параметров. Перечислим их здесь: 1. Сила тока I, протекающего через ТЭОБ; 2. Граничный поток тепла на «горячей» поверхности ТЭОБ (числовые значения потока Qsrm и QSft и форма его профиля: параболическая или «псевдолинейная»); 3.
Толщина диэлектрического слоя AJiel (между алюминиевыми основаниями и припоем); 4. Теплопроводность fC среды в зазорах ТЭОБ и во внешнем пространстве (воздух или вакуум); 5. Ширина Атог зазоров в ТЭОБ. Приведем в сводной таблице изменяющиеся параметры модельной задачи, при которых проводились наши расчеты. В качестве примера результатов расчетов рассмотрим модельную задачу, в которой температура охлаждающей поверхности изменяется по параболическому типу от 290 К на краях до 260 К в середине термоохлаждающей батареи. Профиль потока тепла, отводящегося от горячей поверхности, изменяется по параболе от 3 Вт/см на левом краю до 2 Вт/см в середине; затем поток не изменяется вплоть до правого края, оставаясь при этом равным 2 Вт/см . Поясним, что распределение температуры по горизонтали для наглядности иллюстрируется в пяти различных сечениях (рис.17): Сечение 1: нижняя граница полупроводниковых N- и Р-столбцов с припоем (полудой). Сечение 2: внутренняя часть полупроводниковых N- и Р-столбцов. Сечение 3: верхняя граница полупроводниковых N- и Р-столбцов с припоем (полудой). Сечение 4: Верхняя граница шины - контакта с диэлектрическим слоем, отделяющим алюминиевую подложку. Сечение 5: Верхняя часть алюминиевой подложки у «горячей» поверхности термоохладителя.
Анализ характеристик производительности термоохлаждающей батареи
В работе предложена методика решения уравнения теплового баланса и исследуется двумерное распределение температур в термоэлектрической охладительной батарее, состоящей из большого количества последовательно соединенных термопар. Распределение температуры рассчитывается в области, состоящей из множества подобластей, отличающихся различными физическими свойствами и нелинейной зависимостью коэффициентов, входящих в уравнение теплового баланса.
Расчет всей термоохладительной батареи как единого объекта требует большого объема вычислительных ресурсов и памяти ЭВМ ввиду ее значительной длины, что в большинстве случаев вызывает серьезные трудности в вопросах моделирования. Поскольку описываемая батарея состоит из большого количества одинаковых элементов, то настоящая работа решает задачу одновременного расчета этих элементов на многопроцессорной ЭВМ.
В настоящей работе предложен итерационный метод разделения подобластей со стыковкой на границах решения, учитывающий специфику и особенности задачи. Главный итог исследования - разработка вычислительной методики, позволяющей моделировать тепловые процессы в термоэлектрической батарее с использованием многопроцессорного компьютера с распараллеливанием вычислений в различных отдельных элементах.
Основное решаемое уравнение здесь — двумерное уравнение теплового баланса в различных частях ТЭОБ, которое состоит из нескольких слагаемых. Эти слагаемые учитывают тепловую конвекцию вследствие теплопроводности среды, выделение джоулева тепла в токопроводящих участвках, а также поглощение тепла в полупроводниковых ветвях в результате термоэлектрических явлений Пельтье и Томсона. «Тепло Томсона» имеет характер распределенного в полупроводнике стока тепла. Что касается «эффекта Пельтье», то экспериментальные данные свидетельствуют о том, что тепло Пельтье в рассматриваемой задаче вносит значимый вклад в тепловой процесс лишь на границах полупроводниковых элементов с материалом припоя, прилегающего к ним. При этом на спае, располагающемся с «холодной» стороны полупроводникового элемента, тепло Пельтье поглощается, а на спае с «горячей» стороны тепло Пельтье, соответственно, выделяется.
В основе численных расчетов в каждом типе полупроводниковых столбцов лежит метод приближенной факторизации с разделением двумерной задачи на две одномерные (разновидность метода переменных направлений). Оба этих одномерных уравнения решаются численно методом скалярной прогонки. Так как коэффициенты теплопроводности многократно терпят разрыв, то для решения используется специальная разновидность метода - потоковая прогонка в вертикальном направлении. Такое чередование одномерных прогонок повторяется до тех пор, пока сходимость расчетной задачи не даст величину невязки, равную необходимой . Теоретическая сходимость на практике, однако, осложняется тем, что величину шага разностной сетки невозможно снижать до предельно малых размеров, и поэтому практически достижимая относительная невязка достигает порядка 10 -10 .
Области воздуха (вакуума) между полупроводниковыми ветвями, а также прилегающие к ним области припоя и шины — контакта являются допустимыми для аналитического вычисления решения уравнения теплового баланса. Здесь происходит подгонка решений (сшивка) по горизонтальным границам раздела различных сред до получения необходимой точности решения (малости невязки).
Решения в каждой из расчетных областей стыкуются друг с другом, причем если эти области относятся к сфере расчетов разных процессоров, то процессоры передают данные на соседние в ходе расчетов. Передача данных организована таким образом, чтобы процесс численных расчетов сходился к гладкому решению, и для этого подобраны специальные итерационные параметры.
Проведены расчеты двумерных термоэлектрических охлаждающих батарей с различными показателями коэффициентов теплопроводности композиционных материалов, различными геометрическими размерами слоев, различными граничными условиями теплоотвода и силой электрического тока, протекающего через термобатарею. Исследована сходимость численного метода, скорость расчетов в зависимости от числа процессоров компьютера, задействованных в расчетах одновременно, а также эффективность работы термоохлаждающей батареи в зависимости от интенсивности теплоотвода, т.е. потока тепла па «горячей» поверхности термобатареи.
Исследованы также области наиболее высоких температурных градиентов, могущих влиять на термоупругие нагрузки в батарее.
Для удобства работы пользователей программа снабжена удобным пользовательским интерфейсом, позволяющим легко изменять входные параметры в пользовательском меню и получать нужные распределения искомых функций, линий температуры, графиков невязки - в наглядном графическом виде или на внешнем печатающем приборе. Перспективным выглядит распространение выбранной методики на трехмерную задачу с добавлением элементов модели переноса носителей зарядов (с учетом процессов генерации и рекомбинации электронов и дырок). Таким образом, задача приобретет законченный вид: в качестве исходных условий будет выступать не ток, протекающий через термобатарею, а напряжение на входящих контактах. Практическое внедрение разработанного метода перспективно в термоэлектрических модулях на металлическом основании, выпускаемых, например, Открытым Акционерным Обществом «Импульс».
