Программирование метода конечных элементов.

Программа для решения уравнения Лапласа, использующая формулировку в глобальной системе координат.

Рис. 8 Двумерная теплопередача в брусе квадратного сечения.

Рассматривается двумерная задача о передаче тепла теплопроводностью в случае бруса квадратного сечения, показанного на рис. 8. На двух боковых сторонах бруса задано линейное распределение температуры, а на двух других предполагается идеальная теплоизоляция. Требуется определить линии постоянной температуры (изотермы) внутри бруса, используя метод конечных элементов.

Область ОABCразбивается на 50 трехузловых треугольных элементов с общим числом узлов, равным 36, так, как это

Рис. 9. Разбиение области на 50 конечных элементов.

Таблица 7

Координаты узлов

Узел	Координаты		Узел	Координаты			Узел	Координаты
	x	y		x		y		x	y
1 2 3 4 5 6 7 8 9 10 11 12	0,0 2,0 4,0 6,0 8,0 10,0 0,0 2,0 4,0 6,0 8,0 10,0	0,0 0,0 0,0 0,0 0,0 0,0 2,0 2,0 2,0 2,0 2,0 2,0	13 14 15 16 17 18 19 20 21 22 23 24	0,0 2,0 4,0 6,0 8,0 10,0 0,0 2,0 4,0 6,0 8,0 10,0	4,0 4,0 4,0 4,0 4,0 4,0 6,0 6,0 6,0 6,0 6,0 6,0		25 26 27 28 29 30 31 32 33 34 35 36	0,0 2,0 4,0 6,0 8,0 10,0 0,0 2,0 4,0 6,0 8,0 10,0	8,0 8,0 8,0 8,0 8,0 8,0 10,0 10,0 10,0 10,0 10,0 10,0

Показано на рисунке 9. Основными исходными данными программы являются: а) координаты x и y всех узлов; б) соотношение между локальными и глобальными номерами узлов для всех элементов. Эти данные содержатся в табл 7 и 8 соответственно.

Таблица 8.

Соотношение между глобальными и локальными номерами узлов.

Элемент	Номер узла			Элемент	Номер узла			Элемент	Номер узла
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17	1 1 2 2 3 3 4 4 5 5 7 7 8 8 9 9 10	8 2 9 3 10 4 11 5 12 6 14 8 15 9 16 10 17	7 8 8 9 9 10 10 11 11 12 13 14 14 15 15 16 16	18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34	10 11 11 13 13 14 14 15 15 16 16 17 17 19 19 20 20	11 18 12 20 14 21 15 22 16 23 17 24 18 26 20 27 21	17 17 18 19 20 20 21 21 22 22 23 23 24 25 26 26 27	35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50	21 21 22 22 23 23 25 25 26 26 27 27 28 28 29 29	28 22 29 23 30 24 32 26 33 27 34 28 35 29 36 30	27 28 28 29 29 30 31 32 32 33 33 34 34 35 35 36

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

Так как в этой программе объединение элементных матричных уравнений проводится по узлам, необходимо знать элементы, окружающие каждый узел. Эта информация содержится в массиве NSUR(I, J). Переменная I определяет номер строки и номер рассматриваемого узла одновременно. В столбце J=1 записано число элементов, окружающих узел I, а в столбцах J=1,3,… указаны номера этих элементов . Массив формируется следующим образом. Сначала в первый столбец записываются нули, как это показано в таб. 3.3.

Рис. 10. Блок-схема программы.

Для простоты эта и последующие таблицы содержат данные только для первых десяти узлов. Первым обрабатывается элемент 1, у которого как это видно из табл. 8, глобальными номерами узлов являются числа 1,8 и 7. В массиве общее число элементов, окружающих каждый из указанных узлов, при этом увеличивается на единицу. Номер элемента, а именно1, записывается в следующем столбце J=2 для каждого из этих узлов (табл. 10).

Таблица 9.

Начальное состояние массива NSUR

Узел I	NSUR (I,J)
	Число соседних элементов J=1	Номера соседних элементов
		J=2 J=3 J=4 J=5 J=6 J=7
1 2 3 4 5 6 7 8 9 10	0 0 0 0 0 0 0 0 0 0

Таблица 10.

Состояние массива NSURпосле обработки элемента 1

Узел I	NSUR (I,J)
	Число соседних элементов J=1	Номера соседних элементов
		J=2	J=3 J=4 J=5 J=6 J=7
1 2 3 4 5 6 7 8 9 10	1 0 0 0 0 0 1 1 0 0	1 1 1

Следующий элемент (2) с номерами узлов 1,2 и 8 обрабатывается аналогичным образом. Общее число окружающих узел элементов получается для этих узлов увеличением на единицу, а номер элемента, а именно 2,записывается в следующем столбце, как это иллюстрирует табл. 11. После

Таблица 11

Состояние массива NSURпосле обработки элементов 1 и 2

Узел I	NSUR (I,J)
	Число соседних элементов J=1	Номера соседних элементов
		J=2	J=3	J=4 J=5 J=6 J=7
1 2 3 4 5 6 7 8 9 10	2 1 0 0 0 0 1 2 0 0	1 2 1 1	2 2

Обработки первых 18 элементов указанным способом массив NSUR будет иметь вид, показанный в табл. 12. В результате последовательной обработки всех элементов массива NSUR формируется полностью. Таблица 12

Состояние массива NSURпосле обработки элементов 1-18

Узел I	NSUR (I,J)
	Число соседних элементов J=1	Номера соседних элементов
		J=2	J=3	J=4	J=5	J=6	J=7
1 2 3 4 5 6 7 8 9 10	2 3 3 3 3 1 3 6 6 6	1 2 4 6 8 10 1 1 3 5	2 3 5 7 9 11 2 4 6	4 6 8 10 12 3 5 7	12 14 16	13 15 7	14 16 18

Результаты решения вычерчены в виде изотерм на рис. 11.

Рис. 11. Изотермы в поперечном сечении бруса.

Литература по лекциям

8.1.1. Годунов С.К., Рябенький В.С. Разностные схемы (введение в теорию).

8.1.2. Вазов В., Форсайт Дж. Разностные методы решения уравнений в частных производных.

8.1.3. Гельфорд А.О. Исчисление конечных разностей.

8.1.4. Норри Д., де Фриз Ж. Введение в метод конечных элементов.