Исследование теплового режима в светодиоде с учетом изменения применяемого кристалла

Размещено на http://www.allbest.ru/

Министерство образования и науки Российской Федерации

Федеральное государственное бюджетное образовательное учреждение

высшего профессионального образования

НАЦИОНАЛЬНЫЙ ИССЛЕДОВАТЕЛЬСКИЙ

ТОМСКИЙ ПОЛИТЕХНИЧЕСКИЙ УНИВЕРСИТЕТ

КУРСОВАЯ РАБОТА

На тему «Исследование теплового режима в светодиоде с учетом изменения применяемого кристалла»

По дисциплине: «Энергетика теплотехнологий»

Направление - 130100 Теплоэнергетика и теплотехника

Выполнил студент гр.5БМ6В: Абдраманов Е.Д.

Проверил: к.т.н,

доцент каф. ТПТ Половников В.Ю.

Томск - 2017



Содержание

Введение

1. Основная информация о светодиоде и его конструкции

2. Метод решения задачи

3. Физическая постановка общега вида

4. Листинг программы по решению одномерной задачи

5. Решение нелинейного двумерного уравнения теплопроводности светодиода

6. Блок схема

7. Результаты моделирования

8. Анализ результатов и выводы

Литература

Приложение. Листинг программы при решении двумерной задачи



Введение

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

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

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



1. Основная информация о светодиоде и его конструкции

Светодиод -- это полупроводниковый прибор, трансформирующий электроток в видимое свечение. У светодиода есть общепринятая аббревиатура - LED (light-emitting diode), что в дословном переводе означает "светоизлучающий диод". Светодиод состоит из полупроводникового кристалла (чип) на подложке, корпуса с контактными выводами и оптической системы. Непосредственно излучение света происходит от этого кристала, а цвет видимого излучения зависит от его материала и различных добавок. Как правило, в корпусе светодиода находится один кристалл, но при необходимости повышения мощности светодиода или для излучения разных цветов возможна установка нескольких кристаллов.

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

Светодиод - низковольтный прибор. Обычный светодиод, применяемый для индикации, потребляет от 2 до 4 В постоянного напряжения при токе до 50 мА. Светодиод, который используется для освещения, потребляет такое же напряжение, но ток выше - от нескольких сотен мА до 1А в проекте. В светодиодном модуле отдельные светодиоды могут быть включены последовательно, и суммарное напряжение оказывается более высоким (обычно 12 или 24 В). При подключении светодиода необходимо соблюдать полярность, иначе прибор может выйти из строя. Напряжение пробоя указывается изготовителем и обычно составляет более 5В для одного светодиода.

Яркость светодиода характеризуется световым потоком и осевой силой света, а также диаграммой направленности. Существующие светодиоды разных конструкций излучают в телесном угле от 4 до 140 градусов. Цвет, как обычно, определяется координатами цветности и цветовой температурой, а также длиной волны излучения.

Как реагирует светодиод на повышение температуры? Говоря о температуре светодиода, необходимо различать температуру на поверхности кристалла и в области p-n-перехода. От первой зависит срок службы, от второй - световой выход. В целом с повышением температуры p-n-перехода яркость светодиода падает, потому что уменьшается внутренний квантовый выход из-за влияния колебаний кристаллической решетки. Поэтому так важен хороший теплоотвод. Падение яркости с повышением температуры не одинаково у светодиодов разных цветов. Оно больше у AlGalnP- и AeGaAs-светодиодов, то есть у красных и желтых, и меньше у InGaN, то есть у зеленых, синих и белых[4].

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

Необходимость повышения мощности для увеличения светового потока привела к тому, что традиционная форма корпусного светодиода перестала удовлетворять производителей из-за недостаточного теплоотвода. Надо было максимально приблизить чип к теплопроводящей поверхности. В связи с этим на смену традиционной технологии и несколько более совершенной SMD-технологии (surface montage details - поверхностный монтаж деталей) приходит наиболее передовая технология СОВ (chip on board). Светодиод, изготовленный по технологии СОВ, схематически изображен на рисунке 1.

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

Температурная зависимость

Продолжительная и стабильная работа излучающего диода во многом зависит от эффективного отвода тепла от кристалла. В связи с этим у мощных светодиодов должно быть низкое тепловое сопротивление перехода кристалл-подложка. Например, SMD 5730 и SMD 3014 имеют всего 4°C/Вт, что является достижением современных технологий.

Также нормируются:

максимальная температура p-n-перехода (температура кристалла), которая для SMD приборов может достигать 130°C; температурный диапазон, при котором допускается эксплуатация; температурный диапазон, при котором можно хранить полупроводниковый прибор; температурно-временной график пайки SMD светодиодов[3].



Рисунок-1. Схематическая представление светодиода.

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

2. Метод решения задачи

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

При использовании МКР для задач теплопроводности объект исследования представляют в виде совокупности узлов. Аппроксимируя (заменяя) частные производные ДУ конечными разностями получают систему линейных алгебраических уравнений для определения относительной концентрации, как локальной характеристики в каждом узле сетки. Полученная система является незамкнутой, для ее замыкания используют разностное представление граничных условий. Число уравнений в СЛАУ равно числу неизвестных значений в узлах элементов, на которых необходимо найти решение исходной системы, прямо пропорциональных числу подобластей. Так как каждый из элементов связан с ограниченным количеством соседних, система линейных алгебраических уравнений имеет разреженный вид, что существенно упрощает ее решение. В результате получают замкнутую систему линейных алгебраических уравнений, которую решают численными методами с помощью ЭВМ [1].

3. Физическая постановка общего вида

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



1- металлическая основа (Al+Cu);

2- светизлучающий кристалл из полупроводников.

Рисунок 2. Трехмерный вид рассматриваемых объектов

Физическая постановка для одномерной задачи

Представим конструкцию светодиода как двухслойная пластина состоящих из металла и полупроводникового кристалла с разными теплофизическими характеристиками. В общей конструции нижняя часть металлической основы является адиабатическим поэтому на границе x=0 тепловой поток равен нулю, на границе х=L кристалл контактирует с внешной средой и добавляется нелинейность - излучение. Во внутренней части крисстала существует внутренние линейные потоки тепла. Геометрия задачи показана на рисунке 3.

Рисунок 3. Геометрия задачи

Математическая постановка задачи будет иметь вид:

(1)

где 1 соответствует левой пластине (1 на рис.3),

2 соответствует правой пластине (2 на рис.3).

Начальные и граничные условия можно записать следующим образом:

(2)

(3)

Эту задачу в полной математической постановке будем решать методом конечных разностей на равномерной сетке. Для этого разобьем пластину по толщине на N-1 равных промежутков, т.е. построим конечно-разностную сетку

Рисунок 4. Шаблон конечно разностной сетки

x₂, x₃ ……x_N-1 - координаты внутренних узлов;

x₁, x_N - координаты граничных узлов;

x_i - координата в точке контакта двух сред.

Определим значение температуры в і-ом узле в момент времени t = t_n= n*ф как T (xi ,tn) = T ⁿi . Здесь - шаг интегрирования по временной координат, n - номер шага по времени. Далее заменим дифференциальные операторы в (3) на их конечно- разностные аналоги. Будем пользоваться неявной схемой.

В результате аппроксимации частных производных соответствующими конечными разностями получаем следующую систему линейных алгебраических уравнений (СЛАУ):

Выбранную схему аппроксимации частных производных можно графически представить следующим образом:

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

Полученную систему можно свести к наиболее общему виду:

(5)

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

Предположим, что существуют такие наборы чисел бi и вi (i =1, N =1), при которых

(6)

т.е. трехточечное уравнение второго порядка (5 преобразуется в двухточечное уравнение первого порядка (6). Уменьшим в связи (6) индекс на единицу и полученное выражение T_i-1ⁿ⁺¹=б_i-1T_iⁿ⁺¹+в_i-1 подставим в данное уравнение (5):

откуда получаем



Последнее равенство имеет вид (7) и будет точно с ним совпадать, если при всех i= 2.3,….N-1 выполняются соотношения

(7)

Для определения б_i и в_і по (7) необходимо знать б₁ и в₁ , которые находятся из левого граничного условаия. К тому же в нашей задаче нужно знать прогоночные коэффициенты в точке контакта двух сред б_i* и в_i* .

Далее по формулам (6) последовательно находятся при условии, что найдено из правого граничного условия.

Левое граничное условие:

Отсюда следует, что (8)

Проночные коэффициенты в точке контакта двух сред определяется следующим образом [1]:



(9)

Из правого граничного условия температура на правой границе [1]:

(10)

Исходные данные для вычислении:

- толщину пластины принимаем L= 10 см= 0.1 м

- коэффициент теплопроводности металлической основы(Al-Cu) то есть для первой части пластины л₁ = 164,5 Вт/(м*К);

- коэффициент теплоемкости с_р1 = 883 Дж/(кг*К);

- плотность металлической основы с₁ = 2800 кг/м³ ;

- коэффициент теплопроводности для первого кристалла(GaAs) л₂ =58 Вт/(м*К); для второго образца кристалла(GaP) л₂ =95 Вт/(м*К);

- коэффициент теплоемкости 1) с_р2 = 330 Дж/(кг*К); 2) с_р2 = 484 Дж/(кг*К);

- плотность кристаллов 1) с₂ = 5310 кг/м³ ; 2) с₂ = 3510 кг/м³

- cтепень черноты для обоих кристталов выбираем =0,8.

- коэффициент теплообмена k =500 Вт/(м²*К);

- температура среды T_е = 25 ⁰C.

- температура p-n перехода в полупроводнике GaAs достигает 110-120 ⁰С, а в GaP 130-140 ⁰C, поэтому внутренние источники тепла задаем для GaAs q = 10⁵ Вт/м³, для GaP q =1,5 * 10⁵ Вт/м³.

- Начальную температуру зададим T₀ = 15 ⁰C.

Cоставляем программу решения задачи на язык Pascal ABC.

4. Листинг программы по решению одномерной задачи

program svetodiodizluchenie1x;

uses crt;

const mf=500; sigma=5.669e-8; eps0=1e-5; lamda1=165;

ro1=2800; c1=880; eps=0.8; kapa=500; Te=80; T0=20;

type

vector=array[1..mf] of real;

var {раздел описания переменных, которые мы будем использовать в программе}

i, j, N,N1,N2 : integer;

T,Tn, alfa, beta : vector;

ai, bi, ci, fi, d : real;

a1,a2,h, tau,q: real;

lamda2,ro2,c2: real;

L, t_end, time : real;

f, g : text;

begin clrscr;

{с клавиатуры вводим все необходимые входные параметры}

Writeln('Введите количество промежутков в первой части пластины, N1'); Readln(N1);

Writeln('Введите количество промежутков во второй части пластины, N2');

Readln(N2);

Writeln('Введите коэффиециент теплопроводности кристалла, lamda2');

Readln(lamda2);

Writeln('Введите коэффициент теплоемкости кристалла, с2');

Readln(с2);

Writeln('Введите плотность материала кристалла,ro')

Readln(ro2);

Writeln('Введите внутренние источники тепла кристалла,q')

Readln(q);

Writeln('Введите окончание по времени, t_end'); Readln(t_end);

Writeln('Введите толщину пластины 1, L'); Readln(L);

{определяем общее число узлов в пластине}

N:=N1+N2+1;

{определяем расчетный шаг сетки по пространственной координате}

h:=L/(N-1);

{определяем расчетный шаг сетки по времени}

tau:=t_end/100.0;

{определяем поле температуры в начальный момент времени}

for i:= 1 to N do T[i]:=T0;

{проводим интегрирование нестационарного уравнения теплопроводности}

time:=0;

while time<t_end do {используем цикл с предусловием}

begin

{увеличиваем переменную времени на шаг ф}

time:=time+tau;

{запоминаем поле температуры на предыдущем временном слое}

for i:=1 to N do

Tn[i]:=T[i];

{определяем alfa начальный прогоночный коэффициент на основе левого граничного условия, используя соотношение (8)}

alfa[1]:=1;

beta[1]:=0;

{цикл с параметром для определения прогоночных коэффициентов по формуле (7) в первой части пластины}

for i:= 2 to N1 do

begin

ai:=lamda1/sqr(h);

bi:=2.0*lamda1/sqr(h)+ro1*c1/tau;

ci:=lamda1/sqr(h);

fi:=-ro1*c1*T[i]/tau;

{alfa[i], beta[i] - прогоночные коэффициенты}

alfa[i]:=ai/(bi-ci*alfa[i-1]);

beta[i]:=(ci*beta[i-1]-fi)/(bi-ci*alfa[i-1]);

end;

{определяем прогоночные коэффициенты на границе раздела двух частей, используем соотношения (9)}

alfa[N1+1]:=2.0*a1*a2*tau*lamda2/(2.0*a1*a2*tau*(lamda2+lamda1*(1-alfa[N1]))+sqr(h)*(a1*lamda2+a2*lamda1));

beta[N1+1]:=(2.0*a1*a2*tau*lamda1*beta[N1]+sqr(h)*(a1*lamda2+a2*lamda1)*T[N1+1])/(2.0*a1*a2*tau*(lamda2+lamda1*(1-alfa[N1]))+sqr(h)*(a1*lamda2+a2*lamda1));

{цикл с параметром для определения прогоночных коэффициентов по формуле (8) во второй части пластины}

for i:= N1+2 to N-1 do

begin

ai:=lamda2/sqr(h);

bi:=2.0*lamda2/sqr(h)+ro2*c2/tau;

ci:=lamda2/sqr(h);

fi:=-ro2*c2*T[i]/tau - q;

{alfa[i], beta[i] - прогоночные коэффициенты}

alfa[i]:=ai/(bi-ci*alfa[i-1]);

beta[i]:=(ci*beta[i-1]-fi)/(bi-ci*alfa[i-1]);

end;

{цикл с постусловием, позволяющий итерационно вычислить значение температуры на правой границе, вследствие наличия нелинейности в этом граничном условии}

repeat d:=T[N];

T[N]:=(ro2*c2*sqr(h)*Tn[N]+2.0*tau*(lamda2*beta[N- 1]+kapa*h*Te+eps*sigma*h*(sqr(sqr(Te))-sqr(sqr(d)))))/(ro2*c2*sqr(h)+2.0*tau*(lamda2*(1-alfa[N-1])+kapa*h));

until abs(d-T[N])<=eps0; {значение температуры справа определили}

{используя соотношение (7) определяем неизвестное поле температуры}

for i:= N-1 downto 1 do

T[i]:=alfa[i]*T[i+1]+beta[i];

end; {цикл с предусловием окончен}

{выводим результат в файл}

Assign(f,'res.txt'); Rewrite(f);

Writeln(f,'Толщина пластины L = ',L:6:4);

Writeln(f,'Толщина первой части пластины = ',N1*h:6:4);

Writeln(f,'Толщина второй части пластины = ',N2*h:6:4);

Writeln(f,'Число промежутков по координатев первой части пластины N1 = ',N1);

Writeln(f,'Число промежутков покоординате вовторой части пластины N2 = ',N2);

Writeln(f,'Общее число узлов по координате N = ',N);

Writeln(f,'коэффициент теплопроводности кристалла, lamda2 ',lamda2:6:4);

Writeln(f,'коэффициент теплоемкости кристалла, с2 ',с2:6:4);

Writeln(f,'плотность материала кристалла, ro2 ',ro2:6:4);

Writeln(f,'внутренние источники тепла кристалла, q ',q:6:4);

Writeln(f,'Результат получен с шагом по координате h = ',h:6:4);

Writeln(f,'Результат получен с шагом по времени tau = ',tau:6:4);

Writeln(f,'Температурное поле в момент времени t = ',t_end:6:4);

close(f);

Assign(g,'tempr.txt'); Rewrite(g);

for i:=1 to N do

writeln(g,' ',h*(i-1):10:8,' ',T[i]:8:5); close(g);

end.

Результаты программной реализации представлены на рисунках 12-13.

5. Решения двумерного нелинейного уравнения теплопроводности

Физическая постановка и допущении. При решении двумерного уравнения теплопроводности будем решать тот слой в котором соприкасается поверхности криссталла и металлической подложки Во многих конструкциях светодиода кристалл распологается несимметрично на поверхности подложки, поэтому при решении задачи рассматриваем слой целиком. Геометрия задачи (рис 6) представляется как неоднородная пластина внурий которой расположен полупроводниковый кристалл. В кристалле есть внутренние постоянные тепловые источники. В качестве основного материала выступает сплав (Al+Cu). Принимаем размеры пластины исходя из типичных конструкции светодиода L=H= 0.1 м. Горизантальные границы являются адиабатическими, на вертикальных границах выполняется граничное условия 3 рода с нелинейным излучением.

Рисунок 6. Область решения

Математическая постановка задачи будет иметь вид:

(11)

Начальные и граничные условия запишутся следующим образом:

Т.е. вся расчетная область (рис. 6) покрывается сеткой (рис. 7).



Рисунок 7. Разностная сетка области решения

Введем следующее обозначение: T(xi, yj, tn)=Ti,jn

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

(12)

(13)

Разностные уравнения (12), (13) сводятся к стандартному трехдиагональному виду и решаются последовательно методом прогонки (пункт 2.1.). Сначала для всей области решается уравнение (12), после того как его решение будет найдено, переходят к решению уравнения (13).

Рассмотрим решение уравнения (32) методом прогонки. Приведем это уравнение к виду Тогда коэффициенты А_і, С_і, В_і примут вид:

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

на промежуточном временном слое по формулам (6). После этого приступают к решению уравнения (13). Этапы решения уравнения (13) аналогичны решению уравнения (12).

Определим начальные прогоночные коэффициенты б₁ и в₁ из соотношения Из левого граничного условия следует:

Введем обозначение , тогда



(14)

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

где s - номер итерации, - точность вычислений.

Из правого граничного условия определяем темперутуру T_N^n+1/2 :



(15)

В результате получили нелинейное уравнение (15) для определения температуры на правой границе. Это уравнение также можно решить наиболее простым методом - методом простых итераций.

Исходные данные для вычислении:

- размеры неоднородной пластины принимаем L= H= 10 см= 0.1 м

- коэффициент теплопроводности металлической основы(Al-Cu) то есть для основного материала пластины л₁ = 164,5 Вт/(м*К);

- коэффициент теплоемкости с_р1 = 883 Дж/(кг*К);

- плотность металлической основы с₁ = 2800 кг/м³ ;

- коэффициент теплопроводности для первого кристалла(GaAs) л₂ =58 Вт/(м*К); для второго образца кристалла(GaP) л₂ =95 Вт/(м*К);

- коэффициент теплоемкости 1) с_р2 = 330 Дж/(кг*К); 2) с_р2 = 484 Дж/(кг*К);

- плотность кристаллов 1) с₂ = 5310 кг/м³ ; 2) с₂ = 3510 кг/м³

- cтепень черноты для металлической основы выбираем ₁=₂=0,9.

- коэффициент теплообмена k₁=k₂=500 Вт/(м²*К);

- температура среды T_e1 =T_e2 = 25 ⁰C.

- температура p-n перехода в полупроводнике GaAs достигает 110-120 ⁰С, а в GaP 130-140 ⁰C, поэтому внутренние источники тепла задаем для GaAs q = 10⁵ Вт/м³, для GaP q =1,5 * 10⁵ Вт/м³.

- Начальную температуру области решения зададим T₀ = 15 ⁰C

Составляем программу решения задачи на языке Pascal ABC(приложения).

6. Блок схема



7. Результаты моделирования

В текстовом файле `resulsvet2x-izluch.txt' записывался все входные параметры:

- Результат получен с шагом по координате x hx = 0.01010

- Результат получен с шагом по координате y hy = 0.01010

- Результат получен с шагом по времени tau = 1.8000; 6.00,

- Температурное поле в момент времени t = 180.00 ; 600.00

Текстовый файл `temprsvet2x-izluch.txt' - файл с численными значениями температур в пластине с заданными шагами по двум координатам. Распределение температуры было визуализировано с помощью программы OriginPro и представлено на рисунках 8-11.

При применении кристалла из полупроводников - GaAs(Арсенид Галлий):

Рисунок 8. Температурное поле в момент времени t= 180 c



Рисунок 9. Температурное поле в момент времени t=600 c

светодиод нелинейный одномерный

При применении кристалла из полупроводников - GaP(Фосфид Галлий):

Рисунок 10. Температурное поле в момент времени t=180 c



Рисунок 11. Температурное поле в момент времени t =600 c

Результаты моделирования одномерной задачи:

Рисунок 12. Распределения температуры по толщине металлической основы и кристалла с применением полупроводников GaAs(Арсенид Галлий)



Рисунок 13. Распределения температуры по толщине металлической основы и кристалла с применением полупроводников GaP (Фосфид Галлий).

8. Анализ результатов и выводы

Исследовано нестационарный тепловый режим в светодиоде. В результате выполнения курсового проекта была разработана программная реализация математической модели теплопереноса в неоднородном светодиоде на языке Pascal ABC при заданных условиях. Был изучен конечно-разностный метод Самарского А.А. Для решения дифференциальных уравнений был использован метод прогонки. Уравнение теплопроводности было решено в одномерном и в двумерном виде. При решении двумерного дифференциального уравения теплопроводности для аппроксимации был использован локально-одномерный метод Самарского. Рассматривались две разные светоизлучающие кристаллы, первый- GaAs(Арсенид Галлий), второй - GaP(Фосфид Галлия) с разными теплофизическими характеристиками.

По полученным результатам в решении одномерного дифференциального уравения теплопроводности по толщине светодиода (рисунок 12,13) с помощью температурных градиентов заметим что у светодиода с кристаллом GaP нагретый кристалл быстрее передает тепло к металлической основе это связано с тем, что у полупроводника GaP коэффициент теплопроводности 1,5 раза больше чем GaAs.

По полученным результатам в решении двумерного дифференциального уравнения теплопроводности с представлением светодиода как неоднородная пластина (рисунок 8-11) видим что температурные поля имеет разные изображения. Здесь тоже наблюдается скорость переноса тепла из кристалла GaP в металлическую основу превосходит кристалла GaAs. Интенсивность теплообмена больше у GaP чем у GaAs. Есть еще один недостаток кристалла из полупроводников GaAs: нагреваются на 15-20 ⁰С больше. Поэтому во многих светодиодах применяется полупроводники фосфид галлия GaP.



Литература

1. Кузнецов Г.В. Шеремет М.А. Разностные методы решения задач теплопроводности. - Томск 2007. - 172 с.

2. Самарский А.А. Введение в численные методы: учебное пособие для вузов. Московский государственный университет им. М. В. Ломоносова. -- 3-е изд. СПб.: Лань, 2005. 288 с

3. Интернет сайт: http://duray.ru/database/stati/svetodiod-ustroystvo-printsip-raboty-preimushchestva/

4. Интернет сайт: http://ledjournal.info/spravochnik/ustrojstvo-i-princip-raboty-svetodioda.html

5. Вержбицкий В.М. Основы численных методов. - М.: Высшая школа, 2002. - 840 c.

6. Самарский А.А., Вабищевич П.Н. Вычислительная теплопередача. Москва: Едиториал УРСС, 2003. 784 с.



Приложение. Листинг програмы при решении двумерной задачи

uses crt;

const mf=101; sigma=5.669e-8; eps=1e-5; eps1=0.8;

type

vector1=array[1..mf] of real;

vector2=array[1..mf,1..mf] of real;

var {раздел описания переменных, которые мы будем использовать в программе}

i, j, Nx, Ny : integer;

nx1, nx2 : integer;

ny1, ny2 : integer;

T,Tn: vector2;

lamda1,lamda2,ro1,ro2,c1,c2:real;

a1, a2,d,d1,kapa1,kapa2,Te1,Te2,T0,q : real;

hx, hy, tau, t_end, time : real;

L, H : real;

f, g : text;

procedure progonx(j:integer; lamda,ro,c:real; var W:vector2);

{процедура, разрешающая СЛАУ с трехдиагональной матрицей методом прогонки, в направлении оси Ох ,j - номер слоя по оси у вдоль которого происходит решение СЛАУ;

lamda - коэффициент теплопроводности; ro - плотность; c - коэффициент теплоемкости; W - двумерное поле температуры}

var {раздел описания локальных переменных}

i : integer;

alfa, beta : vector1;

a,ai, bi, ci, fi : real;

begin

a:= lamda/(ro*c);

{определяем alfa начальный прогоночный коэффициент на основе левого граничного условия, используя соотношение (48)}

alfa[1]:=2.0*tau*lamda/(2.0*tau*(lamda+kapa1*hx)+ro*c*sqr(hx));

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

repeat

{определяем beta начальный прогоночный коэффициент на основе левого граничного условия, используя соотношение (48),

при этом начинаем итерационный цикл по левому граничному условию}

d:=W[1,j];

beta[1]:=(ro*c*sqr(hx)*Tn[1,j]+2.0*tau*kapa1*hx*Te1+2.0*tau*eps1*sigma*hx*(sqr(sqr(Te1))-sqr(sqr(d))))/(2.0*tau*(lamda+kapa1*hx)+ro*c*sqr(hx));

until abs(d-W[1,j])<=eps; {значение температуры на левой границе определили}

{цикл с параметром для определения прогоночных коэффициентов по формуле (8)}

for i:= 2 to Nx-1 do

begin

{ai, bi, ci, fi - коэффициенты канонического представления СЛАУ с трехдиагональной матрицей}

ai:=lamda/sqr(hx);

bi:=2.0*lamda/sqr(hx)+ro*c/tau;

ci:=lamda/sqr(hx);

fi:=-ro*c*Tn[i,j]/tau;

{alfa[i], beta[i] - прогоночные коэффициенты}

alfa[i]:=ai/(bi-ci*alfa[i-1]);

beta[i]:=(ci*beta[i-1]-fi)/(bi-ci*alfa[i-1]);

end;

{цикл с постусловием, позволяющий итерационно вычислить значение температуры на правой границе, вследствие наличия нелинейности в этом граничном условии}

repeat

d1:=W[Nx,j];

{определяем значение температуры на правой границе на основе правого граничного условия, используя соотношение (49)}

W[Nx,j]:=(ro*c*sqr(hx)*W[Nx,j]+2.0*tau*(lamda*beta[Nx-1]+kapa2*hx*Te2+eps1*sigma*hx*(sqr(sqr(Te2))-sqr(sqr(d1)))))/(ro*c*sqr(hx)+2.0*tau*(lamda*(1-alfa[Nx-1])+kapa2*hx));

until abs(d1-W[Nx,j])<=eps; {значение температуры на правой границе определили}

{используя соотношение (7) определяем неизвестное поле температуры на промежуточном (n+1/2) временном слое}

for i:= Nx-1 downto 1 do

W[i,j]:=alfa[i]*W[i+1,j]+beta[i];

end; {окончание процедуры progonx}

procedure progony(i:integer; lamda,ro,c:real; var W:vector2);

{процедура, разрешающая СЛАУ с трехдиагональной матрицей методом прогонки, в направлении оси Оу

i - номер слоя по оси x вдоль которого происходит решение СЛАУ;

lamda - коэффициент теплопроводности; ro - плотность; c - коэффициент теплоемкости; W - двумерное поле температуры}

var {раздел описания локальных переменных}

j : integer;

alfa, beta : vector1;

a, ai, bi, ci, fi : real;

begin

{определяется коэффициент температуропроводности}

a:=lamda/(ro*c);

{определяются начальные прогоночные коэффициенты на основе нижнего граничного условия, в соответствии с рассматриваемой постановкой используется соотношение (20)}

alfa[1]:=2.0*a*tau/(2.0*a*tau+sqr(hy));

beta[1]:=sqr(hy)*W[i,1]/(2.0*a*tau+sqr(hy));

{цикл с параметром для определения прогоночных коэффициентов по формуле (8)}

for j:= 2 to Ny-1 do

begin

{ai, bi, ci, fi - коэффициенты канонического представления СЛАУ с трехдиагональной матрицей}

ai:=lamda/sqr(hy); bi:=2.0*lamda/sqr(hy)+ro*c/tau;

ci:=lamda/sqr(hy); fi:=-ro*c*W[i,j]/tau;

{alfa[j], beta[j] - прогоночные коэффициенты}

alfa[j]:=ai/(bi-ci*alfa[j-1]);

beta[j]:=(ci*beta[j-1]-fi)/(bi-ci*alfa[j-1]);

end;

{определяем значение температуры на верхней границе на основе верхнего граничного условия, в нашем случае используется соотношение (21)}

W[i,Ny]:=(lamda*beta[Ny-1])/(lamda*(1-alfa[Ny-1]));

{используя соотношение (7) определяем неизвестное поле температуры}

for j:= Ny-1 downto 1 do

W[i,j]:=alfa[j]*W[i,j+1]+beta[j];

end; {окончание процедуры progonу}

procedure progonxIV(j,n1,n2:integer; lamda,lamdan,ro,ron, c,cn:real; var W:vector2);

{процедура,разрешающая СЛАУ с трехдиагональной матрицей методом прогонки,в направлении оси Оx при наличии включения с отличающимися теплофизическими характеристиками

j - номер слоя по оси у вдоль которого происходит решение СЛАУ; n1 - номер узла с которого начинается включение; n2 - номер узла которым заканчивается включение

lamda - коэффициент теплопроводности основного материала; lamdan - коэффициент теплопроводности материала включения;

ro - плотность основного материала; ron - плотность материала включения; c - коэффициент теплоемкости основного материала;

cn - коэффициент теплоемкости материала включения; W - двумерное поле температуры}

var {раздел описания локальных переменных}

i : integer;

alfa, beta : vector1;

ai, bi, ci, fi : real;

a, an : real;

begin

{определяются коэффициенты температуропроводности}

a:=lamda/(ro*c);

an:=lamdan/(ron*cn);

{определяем alfa начальный прогоночный коэффициент на основе левого граничного условия, используя соотношение (48)}

alfa[1]:=2.0*tau*lamda/(2.0*tau*(lamda+kapa1*hx)+ro*c*sqr(hx));

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

repeat

{определяем beta начальный прогоночный коэффициент на основе левого граничного условия, используя соотношение (48),

при этом начинаем итерационный цикл по левому граничному условию}

d:=W[1,j];

beta[1]:=(ro*c*sqr(hx)*Tn[1,j]+2.0*tau*kapa1*hx*Te1+2.0*tau*eps1*sigma*hx*(sqr(sqr(Te1))-sqr(sqr(d))))/(2.0*tau*(lamda+kapa1*hx)+ro*c*sqr(hx));

until abs(d-W[1,j])<=eps; {значение температуры на левой границе определили}

{цикл с параметром для определения прогоночных коэффициентов по формуле (8) до включения}

for i:= 2 to n1 do

begin

{ai, bi, ci, fi - коэффициенты канонического представления СЛАУ с трехдиагональной матрицей}

ai:=lamda/sqr(hx);

bi:=2.0*lamda/sqr(hx)+ro*c/tau;

ci:=lamda/sqr(hx);

fi:=-ro*c*W[i,j]/tau;

{alfa[i], beta[i] - прогоночные коэффициенты}

alfa[i]:=ai/(bi-ci*alfa[i-1]);

beta[i]:=(ci*beta[i-1]-fi)/(bi-ci*alfa[i-1]);

end;

{определяются прогоночные коэффициенты на границе основного материала и включения}

alfa[n1+1]:=2.0*a*an*tau*lamdan/(2.0*a*an*tau*(lamdan+lamda*(1-alfa[n1]))+sqr(hx)*(a*lamdan+an*lamda));

beta[n1+1]:=(2.0*a*an*tau*lamda*beta[n1]+sqr(hx)*(a*lamdan+an*lamda)*W[n1+1,j])/(2.0*a*an*tau*(lamdan+lamda*(1-alfa[n1]))+sqr(hx)*(a*lamdan+an*lamda));

{цикл с параметром для определения прогоночных коэффициентов по формуле (8) во включении}

for i:= n1+2 to n1+n2 do

begin

{ai, bi, ci, fi - коэффициенты канонического представления СЛАУ с трехдиагональной матрицей}

ai:=lamdan/sqr(hx);

bi:=2.0*lamdan/sqr(hx)+ron*cn/tau;

ci:=lamdan/sqr(hx);

fi:=-ron*cn*W[i,j]/tau-q/2;

{alfa[i], beta[i] - прогоночные коэффициенты}

alfa[i]:=ai/(bi-ci*alfa[i-1]);

beta[i]:=(ci*beta[i-1]-fi)/(bi-ci*alfa[i-1]);

end;

{определяются прогоночные коэффициенты на границе включения и основного материала}

alfa[n1+n2+1]:=2.0*an*a*tau*lamda/(2.0*an*a*tau*(lamda+lamdan*(1-alfa[n1+n2]))+sqr(hx)*(an*lamda+a*lamdan));

beta[n1+n2+1]:=(2.0*an*a*tau*lamdan*beta[n1+n2]+sqr(hx)*(a*lamdan+an*lamda)*W[n1+n2+1,j])/(2.0*an*a*tau*(lamda+lamdan*(1-alfa[n1+n2]))+sqr(hx)*(an*lamda+a*lamdan));

{цикл с параметром для определения прогоночных коэффициентов по формуле (8) в основном материале}

for i:= n1+n2+2 to Nx-1 do

begin

{ai, bi, ci, fi - коэффициенты канонического представления СЛАУ с трехдиагональной матрицей}

ai:=lamda/sqr(hx); bi:=2.0*lamda/sqr(hx)+ro*c/tau;

ci:=lamda/sqr(hx); fi:=-ro*c*W[i,j]/tau-q/2;

{alfa[i], beta[i] - прогоночные коэффициенты}

alfa[i]:=ai/(bi-ci*alfa[i-1]);

beta[i]:=(ci*beta[i-1]-fi)/(bi-ci*alfa[i-1]);

end;

{цикл с постусловием, позволяющий итерационно вычислить значение температуры на правой границе, вследствие наличия нелинейности в этом граничном условии}

repeat

d1:=W[Nx,j];

{определяем значение температуры на правой границе на основе правого граничного условия, используя соотношение (49)}

W[Nx,j]:=(ro*c*sqr(hx)*W[Nx,j]+2.0*tau*(lamda*beta[Nx-1]+kapa2*hx*Te2+eps1*sigma*hx*(sqr(sqr(Te2))-sqr(sqr(d1)))))/(ro*c*sqr(hx)+2.0*tau*(lamda*(1-alfa[Nx-1])+kapa2*hx));

until abs(d1-W[Nx,j])<=eps; {значение температуры на правой границе определили}

{используя соотношение (7) определяем неизвестное поле температуры}

for i:= Nx-1 downto 1 do

W[i,j]:=alfa[i]*W[i+1,j]+beta[i];

end; {окончание процедуры progonxIV}

procedure progonyIV(i,n1,n2:integer; lamda,lamdan,ro,ron,c,cn:real; var W:vector2);

{процедура,разрешающая СЛАУ с трехдиагональной матрицей методом прогонки,в направлении оси Оy при наличии включения с отличающимися теплофизическими характеристиками

i - номер слоя по оси x вдоль которого происходит решение СЛАУ; n1 - номер узла с которого начинается включение; n2 - номер узла которым заканчивается включение;

lamda - коэффициент теплопроводности основного материала; lamdan - коэффициент теплопроводности материала включения;

ro - плотность основного материала; ron - плотность материала включения;

c - коэффициент теплоемкости основного материала; cn - коэффициент теплоемкости материала включения; W - двумерное поле температуры}

Var {раздел описания локальных переменных}

j : integer;

alfa, beta : vector1;

ai, bi, ci, fi : real; a, an: real;

begin

{определяются коэффициенты температуропроводности}

a:=lamda/(ro*c);

an:=lamdan/(ron*cn);

{определяются начальные прогоночные коэффициенты на основе нижнего граничного условия, в данном случае используем соотношение (20)}

alfa[1]:=2.0*a*tau/(2.0*a*tau+sqr(hy)); beta[1]:=sqr(hy)*W[i,1]/(2.0*a*tau+sqr(hy));

{цикл с параметром для определения прогоночных коэффициентов по формуле (8) до включения}

for j:= 2 to n1 do

begin

{ai, bi, ci, fi - коэффициенты канонического представления СЛАУ с трехдиагональной матрицей}

ai:=lamda/sqr(hy); bi:=2.0*lamda/sqr(hy)+ro*c/tau;

ci:=lamda/sqr(hy); fi:=-ro*c*W[i,j]/tau;

{alfa[j], beta[j] - прогоночные коэффициенты}

alfa[j]:=ai/(bi-ci*alfa[j-1]);

beta[j]:=(ci*beta[j-1]-fi)/(bi-ci*alfa[j-1]);

end;

{определяются прогоночные коэффициенты на границе основного материала и включения}

alfa[n1+1]:=2.0*a*an*tau*lamdan/(2.0*a*an*tau*(lamdan+lamda*(1-alfa[n1]))+sqr(hy)*(a*lamdan+an*lamda)); beta[n1+1]:=(2.0*a*an*tau*lamda*beta[n1]+sqr(hy)*(a*lamdan+an*lamda)*W[i,n1+1])/(2.0*a*an*tau*(lamdan+lamda*(1-alfa[n1]))+sqr(hy)*(a*lamdan+an*lamda));

{цикл с параметром для определения прогоночных коэффициентов по формуле (8) во включении}

for j:= n1+2 to n1+n2 do

begin

{ai, bi, ci, fi - коэффициенты канонического представления СЛАУ с трехдиагональной матрицей}

ai:=lamdan/sqr(hy);

bi:=2.0*lamdan/sqr(hy)+ron*cn/tau;

ci:=lamdan/sqr(hy);

fi:=-ron*cn*W[i,j]/tau-q/2;

{alfa[j], beta[j] - прогоночные коэффициенты}

alfa[j]:=ai/(bi-ci*alfa[j-1]);

beta[j]:=(ci*beta[j-1]-fi)/(bi-ci*alfa[j-1]);

end;

{определяются прогоночные коэффициенты на границе включения и основного материала}

alfa[n1+n2+1]:=2.0*an*a*tau*lamda/(2.0*an*a*tau*(lamda+lamdan*(1-alfa[n1+n2]))+sqr(hy)*(an*lamda+a*lamdan));

beta[n1+n2+1]:=(2.0*an*a*tau*lamdan*beta[n1+n2]+sqr(hy)*(a*lamdan+an*lamda)*W[i,n1+n2+1])/(2.0*an*a*tau*(lamda+lamdan*(1-alfa[n1+n2]))+sqr(hy)*(an*lamda+a*lamdan));

{цикл с параметром для определения прогоночных коэффициентов по формуле (8) в основном материале}

for j:= n1+n2+2 to n1+n2 do

begin

{ai, bi, ci, fi - коэффициенты канонического представления СЛАУ с трехдиагональной матрицей}

ai:=lamda/sqr(hy);

bi:=2.0*lamda/sqr(hy)+ro*c/tau;

ci:=lamda/sqr(hy);

fi:=-ro*c*W[i,j]/tau;

{alfa[j], beta[j] - прогоночные коэффициенты}

alfa[j]:=ai/(bi-ci*alfa[j-1]);

beta[j]:=(ci*beta[j-1]-fi)/(bi-ci*alfa[j-1]);

end;

{определяем значение температуры на верхней границе на основе верхнего граничного условия, в данном случае используем соотношение (21)}

W[i,Ny]:=(lamda*beta[Ny-1])/(lamda*(1-alfa[Ny-1]));

{используя соотношение (7) определяем неизвестное поле температуры}

for j:= Ny-1 downto 1 do

W[i,j]:=alfa[j]*W[i,j+1]+beta[j];

end; {окончание процедуры progonуIV}

begin

clrscr;

{с клавиатуры вводим все необходимые входные параметры}

Writeln('Введите длину пластины, L'); Readln(L);

Writeln('Введите толщину пластины, H'); Readln(H);

Writeln('Введите количество пространственных узлов в пластине по оси х, Nx'); Readln(Nx);

Writeln('Введите количество промежутков междуграницей х=0 и включением, nx1'); Readln(nx1);

Writeln('Введите количество промежутков во включении по оси х, nx2'); Readln(nx2);

Writeln('Введите количество пространственных узлов в пластине по оси y, Ny'); Readln(Ny);

Writeln('Введите количество промежутков между границей у=0 и включением по оси y, ny1'); Readln(ny1);

Writeln('Введите количество промежутков во включении по оси y, ny2'); Readln(ny2);

Writeln('Введите коэффициент теплопроводности основного материала,lamda1'); Readln(lamda1);

Writeln('Введите коэффициент теплопроводности материала включении,lamda2'); Readln(lamda2);

Writeln('Введите коэффициент теплоемкости основного материала,c1'); Readln(c1);

Writeln('Введите коэффициент теплоемкости материала включении,c2'); Readln(c2);

Writeln('Введите плотность основного материала,ro1'); Readln(ro1);

Writeln('Введите плотность материала включении,ro2'); Readln(ro2);

Writeln('Введите температуру окруажающей среды на границе х=0,Te1'); Readln(Te1);

Writeln('Введите коэффицент теплообмена на границе х=0,kapa1'); Readln(kapa1);

Writeln('Введите температуру окруажающей среды на границе х=L,Te2'); Readln(Te2);

Writeln('Введите коэффицент теплообмена на границе х=L,kapa2'); Readln(kapa2);

Writeln('Введите плотность теплового потока во включении, q'); Readln(q);

Writeln('Введите начальную температуру ,T0'); Readln(T0);

Writeln('Введите окончание по времени, t_end'); Readln(t_end);

{определяем расчетные шаги сетки по пространственным координатам}

hx:=L/(Nx-1); hy:=H/(Ny-1);

{определяем расчетный шаг сетки по времени}

tau:=t_end/100.0;

{записываем все входные параметры}

Assign(f,'resulsvet2x-izluch.txt');

Rewrite(f);

Writeln(f,'Длина пластины L = ',L:6:4);

Writeln(f,'Расстояние между границей х=0 и включением nx1 = ',nx1*hx:6:4);

Writeln(f,'Длина включения nx2 = ',nx2*hx:6:4);

Writeln(f,'Толщина пластины H = ',H:6:4);

Writeln(f,'Расстояние между границей y=0 и включением ny1 = ',ny1*hy:6:4);

Writeln(f,'ширина включения ny2 = ',ny2*hy:6:4);

Writeln(f,'Число узлов по пространственной координате x в пластине Nx = ',Nx);

Writeln(f,'Число промежутков между границей х=0 и включением = ',nx1);

Writeln(f,'Число узлов по пространственной координате y в пластине Ny = ',Ny);

Writeln(f,'Число промежутков между границей y=0 и включением = ',ny1);

Writeln(f,'Число промежутков во включении = ',ny2);

Writeln(f,'Коэффициент теплопроводности материала пластины lamda1= ',lamda1:6:4);

Writeln(f,'Плотность материала пластины ro1 = ',ro1:6:4);

Writeln(f,'Теплоемкость материала пластины с1 = ',c1:6:4);

Writeln(f,'Коэффициент теплопроводности материала включения lamda2 = ',lamda2:6:4);

Writeln(f,'Плотность материала включения ro2 = ',ro2:6:4);

Writeln(f,'Теплоемкость материала включения с2 = ',c2:6:4);

Writeln(f,'Температура окружающей среды на границе х=0 Te1= ',Te1:6:4);

Writeln(f,'Температура окружающей среды на границе х=L Te2= ',Te2:6:4);

Writeln(f,'коэффициент теплообмена на границе х=0 kapa1= ',kapa1:6:4);

Writeln(f,'коэффициент теплообмена на границе х=L kapa2= ',kapa2:6:4);

Writeln(f,'Начальная температура T0 = ',T0:6:4);

Writeln(f,'Результат получен с шагом по координате x hx = ',hx:6:4);

Writeln(f,'Результат получен с шагом по координате y hy = ',hy:6:4);

Writeln(f,'Результат получен с шагом по времени tau = ',tau:6:4);

Writeln(f,'Температурное поле в момент времени t = ',t_end:6:4);

close(f);

{определяем коэффициент температуропроводности}

a1:=lamda1/(ro1*c1); a2:=lamda2/(ro2*c2);

{определяем поле температуры в начальный момент времени}

for i:= 1 to Nx do

for j:= 1 to Ny do

T[i,j]:=T0;

{проводим интегрирование нестационарного уравнения теплопроводности}

time:=0;

while time<t_end do {используем цикл с предусловием}

begin

{увеличиваем переменную времени на шаг ф}

time:=time+tau;

{запоминаем поле температуры на n-ом временном слое}

for i:=1 to Nx do for j:=1 to Ny do Tn[i,j]:=T[i,j];

{ СЛАУ в направлении оси Ох в промежутке между границей y=0 и включением }

for j:=1 to ny1 do

progonx(j,lamda1,ro1,c1,T);

{на нижней границе включения }

progonxIV(ny1+1,nx1,nx2,lamda1,0.5*(lamda1+lamda2),ro1,0.5*(ro1+ro2),c1,0.5*(c1+c2),T);

{во включении }

for j:=ny1+2 to ny1+ny2 do

progonxIV(j,nx1,nx2,lamda1,lamda2,ro1,ro2,c1,c2,T);

{на верхней границе включения }

progonxIV(ny1+ny2+1,nx1,nx2,lamda1,0.5*(lamda1+lamda2),ro1,0.5*(ro1+ro2),c1,0.5*(c1+c2),T);

{между включением и границей y = H}

for j:=ny1+ny2+2 to Ny do

progonx(j,lamda1,ro1,c1,T);

{ СЛАУ в направлении оси Оy в промежутке между границей x = 0 и включением }

for i:=2 to nx1 do

progony(i,lamda1,ro1,c1,T);

{на левой границе включения 2}

progonyIV(nx1+1,ny1,ny2,lamda1,0.5*(lamda1+lamda2),ro1,0.5*(ro1+ro2),c1,0.5*(c1+c2),T);

{во включении }

for i:=nx1+2 to nx1+nx2 do

progonyIV(i,ny1,ny2,lamda1,lamda2,ro1,ro2,c1,c2,T);

{на правой границе включения }

progonyIV(nx1+nx2+1,ny1,ny2,lamda1,0.5*(lamda1+lamda2),ro1,0.5*(ro1+ro2),c1,0.5*(c1+c2),T);

{в промежутке между включениями и границей x= L}

for i:=nx1+nx2+2 to Nx-1 do

progony(i,lamda1,ro1,c1,T);

end; {цикл с предусловием окончен}

{выводим результат в файл}

Assign(g,'temprsvet2x-izluch.txt');

Rewrite(g);

for i:=1 to Nx do

for j:=1 to Ny do

writeln(g,' ',hx*(i-1):10:8,' ',hy*(j-1):10:8,' ',T[i,j]:8:5);

close(g);

Assign(g,'matrsvet2x-izluch.txt');

Rewrite(g);

for i:=1 to Nx do

begin

for j:=1 to Ny do

write (g,' ',T[i,j]:4:2);

writeln (g);

end;

close(g);

end.

Размещено на Allbest.ru

...