Лекция 4-1. МКЭ на билинейных четырехугольных элементах

⇐ ПредыдущаяСтр 4 из 4

Билинейный элемент

Пусть сетка состоит из четырехугольников, см. рис. 1.

Рис. 1 Базисный квадрат и четырехугольный элемент сетки

Введем базисный квадрат и определим преобразование актуальных координат в локальные :

(1)

На базисном квадрате введем базисные функции , и выразим через них преобразование координат :

(2)

Функции удобно строить с помощью интерполянтов Лагранжа

, (**)

а именно

(3)

Дифференцирование

(4)

Интегрирование

(5)

Вычислим элементы матрицы Якоби с помощью формул (2) – (4).

(6)

Теперь можно подсчитать градиент (4) функции , предварительно разложив по базизу в , а именно, записав

. (7)

Сюда входят градиенты базисных функций и вычислим градиенты . Чтобы вычислить их, подставим (**) в (6). Получим

(8)

Итак, (7) + (8) определяют операцию дифференцирования. А подстановка (8) в (2) и затем в (1) позволяет найти все элементы матрицы Якоби. Например,

Точно также найдем

(9)

Тестовая задача Дирихле для уравнения теплопроводности

Будем решать задачу

(10)

Интегральное тождество

(11)

Граничные условия учитываются позже на матричном уровне.

Покроем область сеткой четырехугольных билинейных элементов и на основе (11) построим схему МКЭ, т.е. СЛАУ. Алгоритм точно такой же, как на треугольной сетке. Для простоты будем считать . Заменив интеграл по области суммой интегралов по элементам и разлагая функции и по базису каждого элемента, получим

‑ из первого интеграла слева получится матрица жесткости:

(12)

Распишем скалярное произведение под интегралом:

(13)

Подставим (13) в (12) и получим формулу для коэффициентов матрицы жесткости элемента

(14)

‑ из второго интеграла слева получится матрица масс:

(15)

‑ из интеграла в правой части получится такая же матрица масс:

(16)

После поэлементной сборки получается глобальная система МКЭ

(17)

Глобальные матрицы получаются «суммированием» элементных матриц.

Осталось рассмотреть вопрос о вычислении интегралов в (14) и (15). Точные вычисления, как в случае треугольных элементов, слишком громоздки, поэтому будем применять квадратурные формулы численного интегрирования. Известно много таких формул разной точности. Иногда можно применять простейшую формулу

В большинстве случаев достаточно использовать 6-узловую формулу

(18)

в которой координаты точек интегрирования и веса представлены в таблице.

k
		0.934172358962716	0.491365692888926e0/4
		-0.356822089773090	0.128641208488885e1/4
	0.774596669241483	0.390885162530070	0.761883709085613e0/4
	-0.774596669241483	0.390885162530070	0.761883709085613e0/4
	0.774596669241483	-0.852765377881771	0.349227402025498e0/4
	-0.774596669241483	-0.852765377881771	0.349227402025498e0/4

В заключение на рис. 2 представим результаты тестового расчета для задачи (10) при

‑ точное решение задачи.

Рис. 2. Численное решение, полученное МКЭ на сетке из 5652 билинейных элементов (слева) в сравнении с точным решением задачи (справа)

Вычисление производных от сеточного решения при использовании билинейных элементов

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

Пусть задана ‑ сеточная функция со значениями в узлах МКЭ-сетки 4-угольных билинейных элементов. Требуется найти вектор

(19)

Для решения задачи (19) запишем интегральное тождество

Интегралы по элементам вычислим с помощью базисного квадрата в локальных координатах . Пользуясь формулами (4), (5), получим

(20)

Для вычисления интегралов применяем те же квадратурные формулы, что и при формировании матриц системы МКЭ. После сборки получаем две СЛАУ:

, (21)

где матрицы и получены «суммированием» элементных матриц.

Решение задач (21) дает искомые производные в узлах МКЭ-сетки. Результаты тестового решения показаны на рис. 3. Хотя на карте изолиний заметны погрешности численного решения, величина этой погрешности невелика, что наглядно демонстрирует рис. 4.

Рис. 3. Производная как решение задачи (21) и точное поле (слева)

Рис. 4. Производная в сечении как решение задачи (21) и точная функция

Лекция 5. Введение в М К О

Основы алгоритма

Требуется решить задачу в области , где ‑ дифференциальный оператор, ‑ заданная функция.

Примеры.

1) Стационарное уравнение теплопроводности

(5.1)

2) Нестационарное уравнение теплопроводности. Полудискретизованная неявная схема по времени с шагом

(5.2)

Можно и умножить все на и разделить на .

3) Конвективная теплопроводность, схема явная по конвекции

(5.3)

Как решать? МКО: покрыть непересекающимися конечными объёмами, ‑ конечный объем.

Примеры.

1) 1D.

2) 2D cтруктурированная

3) 2D неструктурированная

Интегрируем уравнение (5.3) по с границей

(5.4)

Здесь ‑ внешняя нормаль к границе .

Консервативность МКО.

Если сложить все уравнение (5.4), то интегралы по всем внутренним границам КО взаимно уничтожатся, т.к. эти границы попарно-смежные, и нормали у пар направлены противоположно. Останутся лишь интегралы по внешним участкам границы . Будем иметь

Окончательно имеем закон сохранения энергии

Вернемся к системе уравнений (5.4), чтобы получить из него СЛАУ относительно узловых значений. Для этого необходима какая-то аппроксимация при вычислении интегралов в (5.4). Обычно, когда с каждым КО связан один узел , полагают

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

⇐ Предыдущая 1 2 34