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

Здесь выражения N_q^[1] представляют ФФ конечного элемента [1] в q–м узле (q=1,2,3). Для их правильного вычисления в формулы (9.11) и ( 9.9) следует подставлять значения глобальных координат. Только в этом случае

выражение (9.24) действительно позволит учесть соответствие индексов элемента [i, j, k] глобальным номерам узлов. Согласно принятой на рисунке (9.7) нумерации узлов Э_q, соответствие [i j k] примет вид:

Э₁: [231], Э₂: [324], Э₃: [534], Э₄: [635], Э₅: [136] (9.25)

Учитывая принятую нумерацию индексов, приходим к следующей совокупности уравнений для всех конечных элементов D – области:

Рис. 9.7

^[1] = N₂^[1]Ф₂ + N₃^[1]Ф₃ + N₁^[1]Ф₁

^[2] = N₃^[2]Ф₃ + N₂^[2]Ф₂ + N₄^[2]Ф₄

^[3] = N₅^[3]Ф₅ + N₃^[3]Ф₃ + N₄^[3]Ф₄

^[4] = N₆^[4]Ф₆ + N₃^[4]Ф₃ + N₅^[4]Ф₅

^[5] = N₁^[5]Ф₁ + N₃^[5]Ф₃ + N₆^[5]Ф₆

Пример 10.1.

Получить ФФ N₆^[4] и N₆^[5] в заданной на рисунке 9.7 D – области.

Решение.

1. В соответствии с выражениями (9.11) и (9.9) запишем общее выражение для ФФ в произвольной точке элемента 4:

   Ni =	Ai+Bix+Ciy	=	{ (XjYk – XkYj + (Yj – Yk) x + (Xk – Xj) y }
   	2A		2A

2. Пользуясь выражением (9.25), запишем соответствие индексов для конечного элемента Э₄ имеем: [i=6j=3k=5]. Заменяя индексы их значениями, получим для N₆^[4] :

   N6[4] =	{ (X3Y5 – X5Y3 + (Y5 – Y3) x + (X3– X5) y }
   	2A

3. Аналогичное соответствие индексов для конечного элемента Э₅ имеем: [i=1j=3k=6]. Заменяя индексы их значениями, получим для N₆^[5]:

N6[5] =	{ (X3Y5 – X5Y3 + (Y3 – Y5) x + (X5– X3) y }
	2A

Последние две формулы показывают, что ФФ N₆^[5] и N₆^[4] – совершенно разные функции, аппроксимирующие заданный функционал соответственно в пределах конечных элементов 5 и 4. Однако, в самом шестом узле обе эти функции принимают единичные значения, поскольку числители обоих формул в этой точке принимают значения определителя (9.8–а), равные 2А.

6. Решение краевых задач методом конечных элементов.

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

6.1. Задача переноса тепла в стержне.

Постановка задачи. Выберем в качестве ИТО одномерный стержень с коэффициентом теплопроводности , показанный на рисунке 11.1-а. Стержень имеет теплоизолированную боковую поверхность. К левому концу стержня подводится тепловой поток заданной интенсивности q (Вт/см²). На правом конце стержня происходит конвективный обмен тепла с коэффициентом теплообмена – h (Вт/см^{2 о}С). Температура окружающей среды – Т_ос (^оС). Поскольку стержень теплоизолирован, потерь тепла через боковую поверхность не происходит. Требуется определить температурное поле вдоль стержня в установившемся режиме.

Известно, что для данной модели распределение температуры внутри стержня описывает следующее дифференциальное уравнение:

		д2T	= 0			(11.1)
		дx2
a)				б)
Рис. 10.1

При этом, поскольку в установившемся режиме в точках приложения (при х=0) и отвода (х=L) тепла тепловая энергия не должна «задерживаться», должны быть соблюдены следующие граничные условия:

• на левом конце стержня (х=0):

  	дT	+ q = 0	(11.2)
  	дx

• на правом конце стержня (х=L):

	дT	+ h (T – TОС) = 0	(11.3)
	дx

Если тепло отводится от стержня, тепловой поток q должен быть положителен, в противном случае – отрицателен.

Исследования методами вариационного исчисления показывают, что с математической точки зрения в интересующем нас установившемся режиме должен достигать минимума следующий функционал:

 =

V



[

дT

]

2

dV +

S

[

QT +

h

(T – TOC)2

]

dS

(11.4)

2

дx

2

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

 =

V



[

дT

]

2

dV +

S1

(qT ) dS +

S2

h

(T – TOC)2 dS

(11.5)

2

дx

2

С физической точки зрения функционал (11.5) моделирует непрерывность теплового потока в установившемся тепловом режиме. Это означает, что в любой момент времени сумма подводимой (через поверхность S₁) к стержню и рассеиваемой им (через поверхность S₂) тепловой энергии равна энергии, сосредоточенной в объеме (V) стержня. В противном случае, не отводимый от стержня избыток тепловой энергии будет продолжать нагревать стержень, что противоречит условию установившегося режима.

Поскольку, с одной стороны, установившийся режим описывается дифференциальным уравнением (11.1) с граничными условиями (11.2 и 11.3), а, с другой стороны, функционал (11.4) достигает минимума именно в установившемся режиме, то минимум функционала (11.4) и является решением ДУ (11.1) с граничными условиями (11.3).

Температура стержня во всех точках сечения S₁ (S₂) одинакова и равна неизвестной пока (но постоянной в стационарном режиме) величине – Т₁ (Т₃). Учитывая, что в данном случае S₁ = S₂ = A и в силу сказанного, выражение (11.5) принимает вид:

qT1	S1	dS +	h	(T3 – TOC)2	S2	dS =	qT1А +	h	(T3 – TOC)2 А	(11.6)
			2					2

Таким образом, исходное уравнение для определения температуры в каждой точке стержня методом МКЭ примет вид:

 =	V		[	дT	]2	dV +	qT1А +	h	(T3 – TOC)2 А	(11.7)
		2		дx				2

Реализация метода МКЭ включает этапы:

1. Определение подобластей (конечных элементов) и их узловых точек. В данном случае, стержень может быть разбит на два одномерных симплекс – элемента, как это показано на рисунке (10.1-б) с узловыми значениями Т₁, Т₂ и Т₃. Температура внутри элементов находится из формул:

T[1] = N1[1] T1 + N2[1] T2 ;

T[2] = N2[2] T2 + N3[2] T3 ;

(11.8)

ФФ здесь согласно (9.5) равны:

N1[1]	=	(X2 – x)	;	N2[1]=	(x – X1)	;
		L[1]			L[1]
N2[2]	=	(X3 – x)	;	N3[2]=	(x – X2)
		L[2]			L[2]

2. Вычисление частных производных, входящих в выражение (11.7):

дT[1]	=	1	(-T1+T2);	дT[2]	=	1	(-T2+T3)	(11.9)
дx		L[1]		дx		L[2]

3. Разделение интеграла в выражении (11.7) на два (по числу подобластей – конечных элементов, выделенных в пункте 1). Необходимость разбиения интеграла продиктована тем, что производная температуры по переменной х (градиент температуры по оси ОХ), входящая под знак интеграла, не является непрерывной в точке Т₃. Учитывая, что dV=Adx, где А – площадь сечения стержня (А₁ = А₂ = А₃ =А), после разделения и подстановки пределов интегрирования получаем выражение:

x2

x3





[

дT

]2

dV =

[1]A[1]



[

дT

]2dx +

[2]A[2]



[

дT

]2dx

(11.10)

2

дx

2

дx

2

дx

V

x1

x2

4., Проведение подстановки (11.9) в (11.10) и интегрирование:

 V		[	дT	]2	dV =	[1]A[1]	[-T1+T2]2 +	[2]A[2]	[-T2+T3]2	(11.11)
	2		дx			2L[1]		2L[2]

5. Выражаем функционал через узловые значения температуры, для чего объединяем выражения (11.7) и (11.11):

 =	C[1]	(-T1+T2)2 +	C[2]	(-T2+T3)2 +qA[1]T1 +	hA[3]	(-T3+TOC)2		(11.12)
	2		2		2

Здесь приняты следующие обозначения: