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

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

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

(1)

Для решения задачи (1) запишем интегральное тождество (в нем ‑ произвольная пробная функция)

(2)

Интегралы по элементам вычислим с помощью базисного треугольника в локальных координатах .

Пусть на задана функция . Ей соответствует функция , причем соответствие определено формулами (3), а обратное – соотношениями (4).

Рис. 1. Базисный треугольник и его отображение на элемент сетки

Матрица Якоби отображения базисного треугольника (см. рис. 1) на произвольный КЭ с вершинами определяется формулами

(3)

Обратное преобразование, очевидно, имеет вид

(4)

где якобиан преобразования равен

. (5)

Дифференцирование произвольной функции

(6)

Интегрирование произвольной функции

. (7)

Пользуясь формулами (3) ‑ (7) и разложением функций вида

в интегральном тождестве (2) по базису

(8)

получим систему уравнений

(9)

Коэффициенты матрицы масс в левой части (9) собираются из элементных (для каждого КЭ ) матриц масс

(10)

Процедура поэлементной сборки глобальной матрицы размерности такова.

1. Все элементы обнуляются.

2. В цикле по элементам для каждого вычисляется элементная матрица (10) и ее элементы добавяются к соответствующим элементам глобальной матрицы. «Соответствующий» определяется соответствием локальных номеров индексов как вершин треугольника и глобальных номеров узлов МКЭ-сетки, которые прописаны в Таблице 2 – Таблице связности, а именно в её строке для элемента . Например, для элемента , см. рис. 1, строка 27 Таблицы связности имеет вид 18, 19, 23; поэтому локальное значение следует добавить к глобальному элементу матрицы масс .

Правая часть глобальной системы (9) – векторы сил ‑ тоже собирается из элементных вкладов ‑ векторов, которые вычисляются на основе следующих соотношений.

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

(11)

на заданный элементный вектор ; индекс указывает на номер КЭ. Заметим, что (11) определяет две матрицы и , соответствующие компонентам вектора . Вычислим матрицу (11) на базисном треугольнике с помощью формул (6), (7).

(12)

При линейном базисе (8) выражения в квадратных скобках не завмсят от координат и выносятся из-под интеграла. Компоненты градиента базисных функций находим дифференцированием:

(13)

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

(14)

Аналогично для матрицы с использованием (13) из (12) получаем

(15)

Поскольку матрицы (14) и (15) имеют специальный вид (у каждой из них все три строки одинаковы), то произведения этих матриц на элементный вектор дадут два вектора с тремя одинаковыми компонентами:

(16)

Здесь индекс указывает на элемент . Из элементных векторов (16) в цикле по всем конечным элементам сетки собираются правые части двух систем уравнений (9).