Формирование матрицы

Метод конечных элементов предлагает следующую схему вычисления матричных коэффициентов W_ij= (w_i , w_j ) , i=1, 2, ... , n; j=1, 2, ... , N.

Область G разбивается на конечные элементы (треугольники), вершинами которых являются точки p_i, i=1, 2, ... , N (рис.1). Следует отметить, что в таком случае граница Г области G аппроксимируется многоугольником с вершинами p_n+1, p_n+2, ... , p_N.

Рис.1. Триангуляция области G

Базисные функции w_i доопределяются линейно в каждом из треугольников, имеющим своей вершиной точку p_i (рис.2).

Рис.2. Линейное доопределение базисной функции w_i

При таком определении базисных функций w_i интеграл

=

будет равен сумме двух интегралов по треугольникам D _{i j k} и D _{i m j} (таким способом мы обозначили треугольники с вершинами p_i, p_j, p_k и p_i, p_m, p_j). Данные треугольники имеют общее основание (p_i, p_j ). Интегралы по другим треугольникам равны нулю. Рис.3 демонстрирует, что действительно, если треугольники не имеют общего основания (p_i,p_j ), то произведение производных от соответствующих базисных функций (а значит и их производных) равно нулю (рис.3а), поскольку области их ненулевых значений не пересекаются, и отлично от нуля для треугольников с общим основанием (p_i,p_j ) (рис.3б).

Рис.3. Базисные функции w_i и w_j для треугольников,

не имеющих общего основания (p_i,p_j ) (а) и

с общим основанием (p_i,p_j ) (б)

Формализация указанных положений заключается в следующем равенстве

= +

+ . (8)

В случае i=j интеграл по области G будет отличен от нуля для всех треугольников с вершиной p_i

= + +

+ + , (9)

поскольку под знаком интеграла стоит произведение функции самой на себя, что в конечном счете означает пересечение областей ненулевых значений для всех треугольников с общей вершиной p_i (см. рис.2).

Вычислим один из интегралов (8), например, по D _{i j k}. Для этого преобразуем подынтегральное выражение:

=( )( )=

=gradw_i gradw_j=const,

где - орты прямоугольной системы координат.

Равенство выражения константе вытекает из линейности базисных функций. Тогда

= grad w_i grad w_j =

= gradw_i gradw_j S_{D ijk},

где S_{D ijk} - площадь треугольника ijk.

Обозначим высоты треугольника ijk, опущенные из вершин p_i и p_j на стороны (p_j, p_k ) и (p_k, p_i ) как h_ijk и h_jki (рис. 4).

Рис.4. Вычисление grad w_i в треугольнике ijk

Вектор grad w_i направлен по высоте h_ijk (что очевидно), а модуль его равен (см. рис.4) |grad w_i|= .

Учитывая определение скалярного произведения для векторов gradw_i и gradw_j, запишем

= S_{D ijk}.

Вычислив аналогичным способом второй интеграл выражения (8), получим

= S_{D ijk}+

+ S_{D imj}. (10)

В случае i=j интеграл по области G будет равен

= S_Dijk+ S_Dikl+ S_Dilm+ S_{Dimj .} (11)

Вычислив таким же образом все остальные матричные элементы системы (7) и решив ее, например, методом Гаусса, определим неизвестные значения потенциалов U(p_j ), j=1, 2, ... , n во внутренних точках области G.

Замечание.Вычисление площадей и высот треугольников, а также углов между ними не представляет большого труда при фиксации координат (x_i,y_i ) точек p_i, i=1, 2, ... , N и использовании соответствующих свойств треугольников.

1.Площадь треугольника ijk может быть вычислена как половина модуля векторного произведения двух сторон (рис.4)

S_Dijk= =0.5|(x_j-x_i )(y_k-y_i )-(y_j-y_i )(x_k-x_i )|.

2.Высота находится из определения площади треугольника как половины произведения основания на высоту

h_ijk= .

2.Угол между высотами равен углу между сторонами, на которые данные высоты опущены, поэтому используя определение скалярного произведения можем определить

cos(h_ijk,h_jki )=cos[ ]=

= .

З А Д А Н И Е

Рассчитать методом конечных элементов (см. формулу (7)) распределение потенциала в плоскости симметрии системы, сечение которой представляет собой прямоугольник шириной 1 и длиной V (V - номер варианта) (рис.5). При этом левая боковая сторона прямоугольника находится под потенциалом 0, правая - под потенциалом 1, а нижняя и верхняя стороны имеют потенциал 0.5. Количество внутренних узлов p_j принять равным 10 и расположить их на линии, параллельной оси х и проходящей через точку с координатой y=0.5. Количество граничных узлов принять равным 24.

Для формирования матричных коэффициентов W_ij использовать следующий алгоритм. Диагональные коэффициенты W_ii рассчитывать по формуле (11). В случае i j коэффициенты W_ij вычислять по формуле (10), если [p_i,p_j] [p₁,p₂] (j=1, 2, ... , n) или [p_i,p_j] [p₁,p_n+1] (j=n+1, n+2, ... , N) и полагать W_ij=0 в противном случае.

Рис. 5. Триангуляция расчетной области

Провести графическое построение зависимости величины потенциала вдоль оси х при значении координаты y=0.5.

СОДЕРЖАНИЕ ОТЧЕТА

1. Краткое изложение теории.

2. Текст программы.

3. Результаты расчетов в цифровом и графическом виде.

КОНТРОЛЬНЫЕ ВОПРОСЫ

1. Получить систему линейных алгебраических уравнений, являющейся базовой для метода конечных элементов.

2. Какие условия накладываются на базисные функции?

3. Показать, что интеграл по области G вырождается в сумму интегралов по треугольникам с общим основанием.

4. Схема вычислений матричных коэффициентов в методе конечных элементов.

ВАРИАНТЫ ЗАДАНИЙ

Варианты отличаются параметром N в условии задания.