Сагдеева Ю.А., Копысов С.П., Новиков А.К. - Ввеление в метод конечных элементов (1061801), страница 3

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

Текст из файла (страница 3)

После интегрирования по частям от функций требуетсятолько непрерывность базисных и весовых функций.В методе Галеркина весовая функция выбирается равной базисной wl =nPNl , ȳ =αm Nm . В первом слагаемом заменять ȳ на сумму не будем (этоm=1слагаемое уйдет позже за счет граничных условий) nP!Z 1dαm NmZ 1nXdNlm=1dx −Nlαm Nm dx =dxdx00m=112=1 Z 1dȳ Nl−Nl dx.dx 00Кусочно-линейные базисные функции удовлетворяют требованию гладкости, так как они непрерывные. Слева выносим коэффициенты αm за знакинтеграла.

ПолучимZ 1 1 Z 1Z 1nXdNl dNmdȳ αmdx −Nl Nm dx = Nl−Nl dx,dx dxdx 0000m=1где l = 1, 2, . . . , n. Вводя обозначенияklmполучим СЛАУu = (α1 , . . . , αn )T = (y1 , . . . , yn )T ,Z 1dNl dNm=− Nl Nm dx, 1 6 l, m 6 n,dx dx01 Z 1dȳ −Nl dx, 1 6 l 6 n,fl = Nldx 00Ku = f,K = [klm ],f = [fl ].Вспомним, что входящие в аппроксимирующие уравнения определенные интегралы могут быть получены простым суммированием их вкладапо каждому элементуZwl RdV =VM ZXe=1wl RdV.VeВклад интеграла по элементу e с узлами i и j можно вычислить в общейформе.

Причем формула для однотипных элементов будет одна и та же.На элементе V e = (xi , xj ) отличными от нуля функциями будут толькофункции Ni , Nj (рис. 2), то есть, если l 6= i, j, то Nl = 0 на (xi , xj ). ОценимMPeвклад произвольного элемента в сумме klm =klm. Получимe=1eklm= 0, если l, m 6= i, j,Z xj dNi dNjeekij = kji =− Ni Nj dx,dx dxxi!2Z xj dNiee2kii = kjj =− (Ni ) dx,dxxi13xj − x,hix − xidNidNj,= −1/hi ,= 1/hi .hidxdx2ZZZ dNi dNj11dNi1dx = −dx = − ,dx = .2(hi ) V ehidxhiVeV e dx dxZZZxj − x x − xihihiNi Nj dx =dx = ,Ni Ni dx = ,hh63eeiiVVNi =Nj =kii = kjj =1hi− ,hi3kij = kji = −1hi− .hi6(4)Элементная матрица для элемента (xi , xj ) имеет вид00eK =00 ·········kiikji···0kijkjj0··· 00  0 = 00 0 000···hi1hi − 3− h1i − h6i···0···− h1i − h6ihi1hi − 300000Вычислив компоненты матрицы элемента простым суммированием по всемэлементам, получим матрицу K.Процесс формирования глобальной матрицы системы K и глобальноговектора правых частей в методе конечных элементов называется ансамблированием (или сборкой).

Матрицу системы принято называть матрицейжесткости.Запишем вид системы, например, для трех элементов и четырех узлов. Предположим, что все элементы имеют равную длину (he = h), тогдаматрица жесткости приобретает вид 1 h − h1 − h600y1h − 3y2 − 1 − h 2( 1 − h ) − 1 − h0h6h3h6, u =  .K= 0y3 − h1 − h6 2( h1 − h3 ) − h1 − h6 1h1hy400−h − 6h − 3Вычислим вклад элемента (xi , xj ) в вектор правых частей fl =MPe=1(отличными от нуля на элементе (xi , xj ) будут вклады при l = i, l = j)xZ xjZ xjxj − x1 (x − xj )2 jhiNi dx =dx = −= ,hihi22xixixi14fleZxjxiNj dx =Zxjxixx − xi1 (x − xi )2 jhidx = = 2.hihi2xiЗаметим, что в точке 0 и 1 не равны нулю только базисные функцииN1 (0) = 1 и Nn (1) = 1.

Элементные векторы правых частей для первогоэлемента, для внутреннего элемента (xi , xj ) (xi , xj 6= a, b), и для последнегоэлемента имеют вид (ненулевые значения стоят в позиции i, j)0 dȳ .. 0− dx (0) − h21 .  h ..− h21− i 1e.2 f =f n−1 = , f = ..hi  ,hn−1..−− 2 2 . hdȳn−1 .. 0(1) −dx20Для примера из трех элементов после сложения всех элементных векторов,получим глобальный вектор вида dȳ− dx (0) − h21 − h1 − h2 22 f = − h2 − h3  .22dȳh3dx (1) − 2Значения производных в первом и последнем элементе вектора правыхчастей неизвестны, но далее вместо первого и последнего уравнений используем уравнения граничных условий y(0) = y1 = 1, y(1) = y4 = 0.

Длясимметричной матрицы системы условия Дирихле следует вносить следующим образом. Отметим, что до внесения граничных условий получаемаяматрица системы вырождена.Учет граничных условий Дирихле с сохранением симметрии матрицы системы.Пусть в МКЭ получена СЛАУ Ku = f с симметричной матрицей Kи необходимо учесть условие Дирихле us = t0 .

Преобразование системыуравнений представляет собой двухшаговую процедуру. Пусть, например,известно значение u5 = t0 ; преобразование сводится тогда к следующему:1. Все коэффициенты пятой строки матрицы, за исключением диагонального, приравниваются к нулю K5j = 0 при j = 1, n и j 6= 5. Диагональный член приравнивается к единице K55 = 1. Пятая компонента f5вектора f заменяется на значение u5 = t0 .2. Из компонент вектора правых частей, кроме пятой компоненты, вычитается произведение пятого столбца матрицы на значение t, затем пя15тый столбец матрицы (кроме диагонального элемента) обнуляется: fj =fj − Kj5 t, Kj5 = 0, j = 1, n, j 6= 5.5Реализация МКЭ в одномерном случаеИспользуемые данные: n — число узлов; m = n − 1 — число элементов; координаты узлов — массив x = (xi ), i = 1, n; элементная матрицаk = (kij ) размера 2 × 2; элементный вектор правых частей b = (bi ), i = 1, n;глобальная матрица жесткости A = (aij ) размера n × n; вектор правых частей f = (fi ), i = 1, n; вектор решения u = (ui ), i = 1, n; c, d — координатыначала и конца отрезка, на котором рассматривается уравнение; us = t —граничное условие Дирихле.Алгоритм1.

Вычисляем координаты узлов:Шаг сетки h = (d − c)/m;Для i = 1, 2, . . . , n x[i] = c + (i − 1)h;2. В цикле по элементам: Для i = 1, 2, . . . , n − 12.1 Формируем элементную матрицу k. (Для примера из предыдущегопараграфа используем формулы (4)).2.2 Проводим процесс сборки — формируем глобальную матрицу жесткости aii = aii + k11 , aii+1 = aii+1 + k12 , ai+1i = ai+1i + k21 , ai+1i+1 =ai+1i+1 + k22 .2.3 Формируем элементный вектор правых частей b.2.4 Проводим процесс сборки — формируем глобальный вектор правыхчастей fi = fi + b1 , fi+1 = fi+1 + b2 .3. Вносим граничное условие Дирихле us = t в матрицу и вектор правыхчастей3.1 ass = 1; fs = t;3.2 Для k = 1, .

. . , n, k 6= s{ ask = 0;fk = fk − aks · t;aks = 0. }4. Решаем систему Au = f .5. Выводим результат.166Линейный треугольный элементЛинейный треугольный элемент представляет собой треугольник с прямолинейными сторонами и тремя узлами i, j, k, по одному в каждой вершине (см. рис. 3). Выведем вид базисных функций треугольного линейного элемента. При работе с узлами требуется определенная нумерацияили направление обхода узлов элемента.

Будем использовать направлениепротив часовой стрелки. Узловые значения в элементе обозначим (ui , uj ,uk ). Координаты узлов соответственно равны (xi , yi ), (xj , yj ), (xk , yk ). Площадь треугольника равна A. Будем интерполировать неизвестную функцию u(x, y) на элементе линейной функцией, то естьue (x, y) = α1 + α2 x + α3 y.(5)Интерполяция проводится так, чтобы значения точного решения и приближенного совпадали в узлах сетки, то есть ui = α1 + α2 xi + α3 yi ,uj = α1 + α2 xj + α3 yj ,uk = α1 + α2 xk + α3 yk .Эта система уравнений всегда имеет единственное решение, так как определитель системы равен двум площадям треугольника, то есть не равеннулю. Разрешив ее относительно αi и подставив αi в (5), получимue =1(ui (xj yk − xk yj + (yj − yk )x + (xk − xj )y) + uj (xk yi − xi yk +2A(yk − yi )x + (xi − xk )y) + uk (xi yj − xj yi + (yi − yj )x + (xj − xi )y).Функции при коэффициентах ui , uj , uk определяют линейные функцииформы данного элемента.

С их помощью соотношение, определяющее элемент, записывается в виде1(al + bl x + cl y)2A aj = xk yi − yk xi  ak = xi yj − xj yibj = yk − yibk = yi − yjcj = xi − xkck = xj − xiue = Ni ui + Nj uj + Nk uk , где Nl = ai = xj yk − xk yjbi = yj − ykci = xk − xjФункция Nl равна единице в узле (xl , yl ) и нулю в двух других узлах.17(6)(xk , yk )kkkhL1L2jiBj(xj , yj )ibL1sL3j(xi , yi )iРис. 3: Треугольный элемент и L-координатыЕсли неизвестная функция аппроксимируется базисными функциями,линейными по x и y, то градиенты в направлениях x и y будут постоянны.Например,∂Ni∂Nj∂Nk∂u=ui +uj +uk ,∂x∂x∂x∂x∂Nlbl=,∂x2Al = i, j, k,X bl∂u=ul = const .∂x2Al=i,j,kЕстественная система координат (ЕСК). При вычислении элементных матриц жесткости или векторов правых частей приходится интегрировать функции формы или их частные производные по площадиэлемента.

Интегрирование можно упростить, если записать интерполяционные соотношения в системе координат, связанной с элементом, то естьв локальной системе координат (ЛСК). Эту локальную систему координат называют естественной, если в узле элемента локальная координатапринимает значения 0, ±1.В двумерном случае интеграл от функции, заданной в глобальной системе координат, может быть вычислен в ЛСК с помощью соотношенияZZf (x, y)dxdy =f (x(r, s), y(r, s)) det |J| dr ds,VV∗где V и V ∗ — соответственно старая и новая области интегрирования,det |J| — абсолютное значение определителя матрицы Якоби преобразования системы координат.

Характеристики

Список файлов книги