Основные положения метода конечных элементов

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

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

К недостаткам метода следует отнести: (а) в основе МКЭ лежит дискретизация всего ИТО, что неизбежно ведет к очень большому количеству конечных элементов (особенно в трехмерных задачах); (б) МКЭ часто приводит к нереальным разрывам значений физических величин между смежными элементами.

При решении задач методом КЭ используются одномерные, 2 и 3-мерные КЭ. Одномерный КЭ показан на рисунке 3.1.1. Площадь поперечного сечения одномерного КЭ может изменяться по длине, но в большинстве практических задач ее считают постоянной. Наиболее часто такой элемент используется в одномерных задачах распространения тепла и в задачах строительной механики.

 

Рис. 3.1.1

Для построения дискретной модели двумерной области используют треугольники и 4-хугольники, которые могут иметь как прямолинейные, так и криволинейные стороны.

Собственно процесс дискретизации ИТО может быть разделен на два этапа: (а) разбиение ИТО на КЭ и (б) нумерация элементов и узлов. Последний этап может существенно повлиять на эффективность вычислений.

Требование простоты КЭ связано с тем, что при моделировании области должно быть использовано большое число элементов, поэтому деление области на треугольники является наилучшим способом разбиения.

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

На практике используются три типа элементов: симплекс-элемент, комплекс-элемент и мультиплекс-элемент.

Симплекс - элементам соответствуют полиномы, содержащие константу и линейные члены. Число коэффициентов в таком полиноме на 1 больше размерности координатного пространства. Например, полином:

j = a₁ +a₂ x +a₃ y

представляет собой симплексную функцию для двумерного треугольного элемента. Этот полином линеен по X и Y и содержит три коэффициента, потому что треугольник имеет три узла.

Комплекс – элементам соответствуют полиномиальные функции, содержащие константу, линейные члены, а также члены второго, третьего и более высоких порядков, если это необходимо. Форма комплекс – элементов может быть такой же как у симплекс - элементов, но с дополнительными граничными (и даже внутренними) узлами. Число узлов в комплекс – элементе должно быть больше размерности координатного пространства + 1. Интерполяционный полином для 2-мерного треугольного комплекс – элемента имеет вид:

j = a₁ + a₂ x + a₃ y + a₄ x ² + a₅ xy + a₆ y ²

Это соотношение включает шесть коэффициентов, поэтому рассматриваемый элемент должен иметь шесть узлов.

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

 

Этапы решения в МКЭ.

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

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]	=		(-T1+T2);	д T[2]	=		(-T2+T3)	(11.9)
дx	L[1]	дx	L[2]

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

x2

x3

ò

l

[

д T

]2

d V =

l[1]A[1]

ò

[

д T

]2dx +

l[2]A[2]

ò

[

д T

]2dx

(11.10)

дx

дx

дx

V

x1

x2

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

ò V	l	[	д T	]2	d V =	l[1]A[1]	[-T1+T2]2 +	l[2]A[2]	[-T2+T3]2	(11.11)
	дx	2L[1]	2L[2]

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

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

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

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

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

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

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

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

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

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

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

C

·

T

=

F

(11.16)

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

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

Решение.

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

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

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

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

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

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

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

	-20			Т1
-20		-20	·	Т2	=
	-20			Т3

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

 

⇐ Предыдущая12345678Следующая ⇒

Поделиться:

Воспользуйтесь поиском по сайту: