Автоматихация производства ЭВА, страница 10

С⁽¹⁾ = (А⁽¹⁾⁽¹⁾/L⁽¹⁾); С⁽²⁾ = (А⁽²⁾⁽²⁾/L⁽²⁾)

6. Получение системы алгебраических уравнений. Правильными значениями Т₁, Т₂ и Т₃ являются те, при которых величина функционала  достигает минимума. Приравнивая нулю первую производную функционала (11.12) по Т₁, получаем первое уравнение системы:

	д	= C[1] T1 - C[1] T2 + qA[1] = 0	(11.13)
	дT1

Аналогично получаем еще два уравнения:

	д	= -C[1] T1 + [C[1] +C[2] ]T2 -C[2] T3 = 0	(11.14)
	дT2
	д	= -C[2] T2 + [C[3] +hA3 ]T3 - hA3TOC = 0
	дT3

Запишем полученную систему в матричной форме:

С(1)	-С(1)	0		Т1		-qA1
-С(1)	С(1)+С(2)	-С(2)		Т2	=	0	(11.15)
0	-С(2)	С(2)+hA3		Т3		hA3TOC

В более общей матричной форме система примет вид:

C



T

=

F

(11.16)

Матрица C в формуле (11.16) называется «глобальной матрицей жесткости ». В контексте задачи переноса тепла –это – «глобальная матрица теплопроводности». Вектор-столбец F называется «глобальным вектором нагрузки ». Искомый вектор [T] будем называть вектором решения.

Пример 11.1. Рассчитать температурное поле в круглом стержне с площадью поперечного сечения A=1 см² и длиной L=7,5 см с теплоизолированными стенками. К левому концу стержня подводится тепловой поток q = 150 Вт/см². Коэффициент теплопроводности материала стержня и коэффициент конвективного теплообмена на правом конце стержня соответственно равны: =75 Вт/(см  ^ОС), h = 10 Вт/(см²  ^ОС). Температура окружающей среды равна Т_ОС=40 ^ОС.

Решение.

1. Тепло подводится к стержню, поэтому тепловой поток q следует записывать со знаком «минус»: q = - 150 Вт/см².

2. Рассчитываем значение термов, входящих в коэффициенты матриц C и F:

С⁽¹⁾ =(А⁽¹⁾⁽¹⁾/L⁽¹⁾)=(175/3,75)=20Вт/(см^ОС),

С⁽²⁾ =(А⁽²⁾⁽²⁾/L⁽²⁾)=(175/3,75)=20Вт/(см^ОС),

hA₃=10Вт/(см^ОС), -qA₁= -(-150)1 = 150Вт/см,

hA₃T_OC=10140 = 400Вт/см.

3. Окончательная система уравнений примет вид:

20	-20	0		Т1		150
-20	40	-20		Т2	=	0
0	-20	30		Т3		400

4. Решением полученной системы являются следующие узловые значения температуры: Т₁=70 ^оС, Т₂=62,5 ^оС ; Т₃=55 ^оС.

Проблема реализации МКЭ на ЭВМ. Процедура минимизации  приводит к системе уравнений, которые решаются относительно узловых значений температур. Однако, с точки зрения реализации процедуры минимизации  на ЭВМ, целесообразно функционал (11.4) представить в виде суммы вида:

m

 = 1 + 2 +…+m =



1

(11.17)

i =1

где: m – количество конечных элементов, на которые разбивается ИТО.

Дело в том, что в библиотеке САПР, реализующей минимизацию функционала  на ЭВМ, содержаться модели не всего ИТО, а именно конечных элементов (например, симплекс – элементов), причем мощность указанной библиотеки КЭ и определяет функциональные возможности САПР ИТО. В процессе решения задачи ЭВМ (в соответствии с заданием на проектирование) автоматически объединяет модели конечных элементов в единую модель ИТО. В этой связи, представляется целесообразным описать последовательность шагов получения системы линейных уравнений (11.16), используя в качестве исходного шага разбиение (11.17). Тем более, что эта процедура и является центральной в работе инженера, моделирующего поведение ИТО на ЭВМ.

Из примера (11.1) ясно, что функционалы по отдельным конечным элементам, выраженные через узловые значения, имеют вид:

1 =





(-T1+T2)2dV

+



qT1 dS

2(L[1])2

V[1]

S[1]

(11.18)

2 =





(-T2+T3)2dV

+



h

(T3+TOC)2dS

2(L[2])2

2

V[2]

S[2]

Проведем дифференцирование (₁) системы (11.18) по всем узловым значениям:

д(1)	=					(-T1+T2) (-1)dV +			q dS
дT1				(L[1])2
		V[1]					S[1]

д(1)	=					(-T1+T2) dV
дT2				(L[1])2
		V[2]

д(1)	= 0
дT3

Вычисляя в этой системе интегралы, и применяя обозначения, принятые в формуле (11.12), получим следующую систему уравнений в обычной и матричной форме:

	д[1]	= + C[1] T1 - C[1] T2 + qA [1]
	дT1
	д[1]	= - C[1] T1 + C[1] T2 + 0
	дT2
	д[1]	= 0 + 0 + 0
	дT3

	д[1]	=	C[1]	-C[1]	0		T1	+	qA[1]
	дT1
	д[1]	=	-C[1]	C[1]	0		T2		0
	дT2
	д[1]	=	0	0	0		T3		0
	дT3

Для краткости изложения будем далее обозначать ее так:

д[1]

д[T]

Запишем систему уравнений (11.19) в матричной форме для первого КЭ:

д(1)	= [ C (1) ]  [ T ] +[ F ]	(11.19)
д[T]

В отличие от системы уравнений (11.16) в системе (11.19) матрица коэффициентов C⁽¹⁾ называется «матрицей жесткости элемента ». Ее название в контексте задачи переноса тепла - «матрица теплопроводности элемента ». Вектор-столбец F как и ранее является «глобальным вектором нагрузки ».

Проведем теперь дифференцирование второй компоненты (₂) системы (11.18) по всем узловым значениям:

д(2)	= 0
дT1

д(2)	=					(-T2+T3)( -1) dV
дT2				(L[2])2
		V[2]

д(2)	=					(-T2+T3) dV +			h (T3-TOC) dS
дT3				(L[2])2
		V[2]					S[2]

После вычисления интегралов получим систему уравнений:

д (2)	=	0	+	0	+	0
дТ1
д (2)	=	0	+	С(2)Т2	-	С(2)Т3
дТ2
д (2)	=	0	-	С(2)Т2	+	(С(2)+hA3)Т3
дТ3

Или в матричной форме:

 (2)	=	0	0	0		Т1		0
 Т1
 (2)	=	0	С(2)	-С(2)		Т2	+	0
 Т2
 (2)	=	0	-С(2)	(С(2)+hA3)		Т3		+hA3
 Т3

Учитывая аддитивный характер функционала  можно утверждать, что для его минимизации по узловым значениям необходимо, чтобы выполнялось равенство:

д	=	д[1]	+	д[2]	= 0	(11.20)
д[T]		д[T]		д[T]

Поэтому, складывая выражения для обеих компонент функционала в матричном виде, получим окончательную систему уравнений:

С(1)	-С(1)	0		Т1		qA1		0
-С(1)	С(1)+С(2)	-С(2)		Т2	+	0	=	0
0	-С(2)	С(2)+hA3		Т3		-hA3TOC		0

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