I_X⁺ = A_X [ ^t_i+1,j,k – ^t_i,j,k] ; I_X^– = A_X [  ^t_i,j,k–  ^t_i-1,j,k] ;

I_Y⁺ = A_Y [ ^t_i,j+1,k – ^t_i,j,k] ; I_Y^– = A_Y [  ^t_i,j,k–  ^t_i,j-1,k] ;

I_Z⁺ = A_Z [ ^t_i,j,k+1 – ^t_i,j,k] ; I_Z^– = A_Z [  ^t_i,j,k–  ^t_i,j,k-1] .

Решение системы уравнений (20) можно находить следующим образом:

1. Определяем температуру в узле через один шаг  по времени:

   t+1i,j,k = I X+ – I X– + I Y+ – I Y– + I Z+ – I Z– + Q +  ti,j,k

   (21)

2. Вычисление по (21) выполняем для всех узлов.

3. Определяем температуру в узлах в момент 2 и так далее.

Если в результате такого расчета разность ^t+1 – ^t будет постепенно уменьшаться и температура стабилизироваться, то решение считается устойчивым. В противном случае – решение не устойчиво.

3.3. Моделирование процесса пайки выводов ЭРЭ

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

Исследования показывают, что воздействие t⁰ при пайке аналогично термоудару на ЭРЭ, последствия которого (отслаивание подложки, нарушение герметичности корпуса и др.) не проявляются сразу при монтаже аппаратуры, а являются причиной ее отказа при эксплуатации. Отсюда следует, что обеспечение нормального теплового режима пайки имеет важное значения для обеспечения надежности всей аппаратуры.

Основными параметрами пайки является температура пайки Т_П, время пайки t_П, температура нагрева прибора Т_Н и расстояние от корпуса прибора до места пайки L.

Нагрев ЭРЭ в результате пайки представляет сложный процесс передачи тепловой энергии от места пайки к корпусу ЭРЭ, в котором участвуют все виды передачи тепла: кондуктивный теплообмен по выводу, излучение и конвективный теплообмен поверхностей ЭРЭ и выводов с окружающей средой. Сам процесс является нестационарным, так как за время пайки происходят изменения температуры в различных точках изделия. Подобные задачи на практике решаются с использованием различных допущений и идеализаций.

Электронный прибор представляют в виде стержня, на одном конце которого расположен источник постоянной температуры (место пайки), а на другом - конструктивно связанная со стержнем фиктивная масса (ФМ). Эта ФМ представляет собой математическую модель реального электронного прибора, обладающую адекватными прибору теплофизическими параметрами: теплоемкостью и теплоотдачей в окружающую среду. Геометрически ФМ можно представить в виде тонкого диска, размещенного на торце внешнего вывода. Теплопроводность ФМ (массивного ЭРЭ) считается бесконечно большой, то есть температура ФМ во всех ее точках одинакова.

Нагрев ФМ аналогичен нагреву ЭРЭ в интересующей точке (например, в месте крепления полупроводникового кристалла в корпусе). Сказанное позволяет свести задачу математического моделирования процесса нагревания ЭРЭ при пайке к задаче определения температуры ФМ, если сделать следующие предположения:

Рис. 6	Рис. 7

• нагрев ФМ равномерен в каждый момент времени;

• конец вывода в месте пайки мгновенно достигаем температуры пайки;

• тепловые коэффициенты не зависят от температуры;

• конвективный обмен вывода с окружающей средой пренебрежимо мал.

При таких допущениях тепловой процесс пайки представим с помощью двух ДУ:

S	дТ	- АТ = СФМ			дТ	(8.1)
	дХ				дt
	д2Т	= C	дТ	(8.2)
	дХ2		дt

где: S – площадь сечения вывода ЭРЭ (S = 10^–6м² –далее в скобках приводятся значения, используемые далее в примере расчета);

 – коэффициент теплопроводности материала вывода (=400 Дж/мс^оК);

– изменение температуры на границе вывода и корпуса ЭРЭ;

А – коэффициент теплоотдачи ЭРЭ во внешнюю среду (А=0,5 Дж/(с^оК) );

С_ФМ – теплоемкость фиктивной массы (С_ФМ =0,05 Дж/^оК);

С– удельная теплоемкость материала вывода (С=410⁶ Дж/(м³^оК);

– скорость изменения температуры ФМ;

X – расстояние от места пайки; T – температура ФМ; t – время пайки.

Нагрев ФМ в различные моменты времени получается решением уравнения (8.1). В нем уменьшаемое и вычитаемое в левой части характеризует подвод тепла к ЭРЭ через вывод и теплоотдачу ЭРЭ в окружающую среду, а правая часть – описывает нагрев ЭРЭ в процессе пайки.

Уравнение (8.2) определяет передачу тепла через вывод ЭРЭ. В нем левая часть задает тепловой поток через вывод, а правая – нагрев вывода.

Начальным условием обоих уравнений является температура вывода и ЭРЭ в нулевой момент времени начала пайки:

T(X,0) = 293^OK (20^OC) (8.3)

Граничным условием на (левом) конце вывода является температура пайки:

T(0,t) = Т_П (8.4)

Для решения уравнений (8.1) и (8.2), относящихся к классу линейных ДУ (1- уравнение первого рода и 2- параболическое уравнение) разделим вывод элемента поперечными сечениями с шагом h (рисунок 7) и на временной оси отметим моменты времени Тj=j где: – интервал времени, j – индекс ( j=0 в начальный момент пайки).

Температура вывода в i–м сечении в j-й интервал времени от начала пайки:

T_ij = T ( ih, j )

Учитывая формулы (9) и (11), полученные в разделе 3, можно записать следующие выражения перехода от частных производных к конечным разностям в соответствии с принятыми обозначениями:

дТij		Тij- Тi-1,j	;	дТij		Тi,j+1 - Тi,j	(8.5)
дХ		h		дt		

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

д2Тij		Тi+1,j- 2Т1j + Тi-1,j	(8.6)
дХ2		h2

Заменим в уравнении (8.1) производные конечными разностями, получим:

S	Тn,j- Тn-1,j	- AТn,j = СФМ	Тn,j+1 - Тn,j	(8.7)
	h		

Отсюда следует, что в момент времени t=(j+1) температура ФМ составит:

Тn,j+1 = Тn,j{	S	-	A	+1} - Тn-1,j	S	(8.8)
	hСФМ		СФМ		hСФМ

Заменим в уравнении (8.2) вторые производные конечными разностями (8.6):

	Тi+1,j–2Тi,j+Тi-1,j	= C	Тi,j+1 – Тi,j
	h2		

Пусть  =(/h²C). Можно показать, что условие устойчивости решения, то есть возможность получения точных решений за определенное количество шагов, примет вид:

	h2C	=	1	 
	2		2

Далее примем =0.2/. Тогда значение рабочей температуры в момент времени t=(j+1) в i –м сечении стержня примет вид:

Ti,j = 0.2Ti+1,j + 0,6 Ti,j + 0.2Ti-1,j

(8.9)

По формуле (8.8) и (8.9) можно рассчитать температуру ФМ в сечениях вывода последовательно в моменты времени t_j = 0, t_j+1 = , t_j+2 = 2, … . В общем случае, в произвольный k-й момент времени t_k = (t_k-1 + ). Программа на языке Паскаль вычисления температурного поля в стержне методом конечных разностей (МКР), действующая по формуле (8.9), может иметь следующий вид:

CONST n=10; Delta=0.1; L=20; Type mass=array[1..n]of LongInt;

Lambda=400; A=0.1; Cfm=6; C=4000000; S=0.000001; h=0.001;

VAR Md,Mt:mass; i,t:integer; b:boolean; k1,k2:Real;

BEGIN

k1:=((S*lambda*Delta)/(h*Cfm))-(delta*A/Cfm)+1; k2:=S*lambda*Delta/(h*Cfm);

For i:=1 to n Do Mt[i]:=293; Mt[1]:=600; For i:=1 to n Do MD[i]:= Mt[i];

Repeat For i:=2 to n-1 Do

Md[i]:=round(0.2*(Mt[i+1]+Mt[i-1])+0.6*Mt[i]);

md[n]:=round(k1*mt[n] - k2*mt[n-1]); For i:=2 to n Do Mt[i]:=Md[i];

Until false; END.

По формуле (8.9) аналогичным образом может быть составлена программа, реализующая вычисление значений температуры в дискретные моменты времени в последнем n-м сечении стержня, к которому подсоединена фиктивная масса.

Ошибки ограничения уменьшаются при h  0 и   0, то есть решение разностных уравнений (8.8) и (8.9) асимптотически приближается или сходится к решению ДУ (8.1) и (8.2). Следовательно, с целью повышения точности решений следует уменьшать шаг дискретизации по длине и по времени, что, в свою очередь, ведет к возрастанию количества основных циклов в алгоритме вычислений (пропорционально m x n), то есть – к увеличению времени счета.

Другим важным свойством МКР является устойчивость решений. Устойчивость метода означает, что любые ошибки (округления или ограничения) не возрастают в ходе вычисления решения. Доказано, что указанный процесс сходится и устойчив, если (<0.5) или (<1/[2]). Тем самым накладываются ограничения на выбор шага  по времени при выбранном шаге h разбиения вывода сечениями.

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

4. Решение системы алгебраических уравнений на ЭВМ

При использовании метода конечных разностей и методоа конечных элементов получается система линейных алгебраических уравнений, которая должна быть решена относительно неизвестных узловых параметров. Методы решения этой системы с малым или большим числом уравнений мало отличаются друг от друга. Специфика ее решения заключается в том, что при дискретизации D-области можно контролировать расположение коэффициентов в глобальной матрице жесткости. В частности, можно получить матрицу ленточного типа, которая имеет два полезных свойства: а) симметричность относительно главной диагонали и б) положительная определенность, означающая, что все коэффициенты на главной диагонали – положительные. Указанные свойства значительно сокращают объем вычислений. Причем минимизируются ошибки округления. Ленточный характер матрицы позволяет значительно сократить объем памяти для ее хранения.

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