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

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

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

Пусть, например, теплообмен задан на стороне13. Вычисляем длину стороны L13 . Тогда к вектору правых частейдобавляем вектор  ZN11N2  dS = hT∞ L13 0 .hT∞2SN31Определение расчетных величин в элементах (пункт 6). Если естьнеобходимость вычислить градиенты температур в узлах элементов  ∂T T11 bi bj bk  ∂xT2 .=∂Tci cj ck2A∂yT326В цикле по элементам1. Определяем координаты узлов текущего элемента и значения решенийДля i = 1 до 3 {xi = mesh[si , 1]; yi = mesh[si , 2]; Ti = u[si ]}.2. Формируем bi , ci :b1 = y2 − y3 ; b2 = y3 − y1 ; b3 = y1 − y2 ;c1 = x3 − x2 ; c2 = x1 − x3 ; c3 = x2 − x1 .1 xi3.

Вычисляем площадь треугольника A = 21 det 1 xj1 xkx3 y2 + x1 y2 − x1 y3 + x3 y1 − x2 y1 )/2.yiyj  = (x2 y3 −yk4. Вычисляем градиенты температур на всем элементе (они будут постоянны по элементу)∂T  T1 ∂x 1b1 b2 b3  T2 .= ∂T  2A c1 c2 c3T3∂yПоскольку градиенты получаются кусочно-постоянными на элементах, можно провести сглаживание значений. Например, если один узел принадлежит четырем элементам, то значению в узле приписывают среднее по этимчетырем элементам.9Применение четырехугольных конечныхэлементов в задаче упругостиПостановка плоской задачи теории упругости. Система разрешающих уравнений теории упругости относительно трех неизвестных полейперемещений u = (ux , uy )T , деформаций ε = (εxx , εyy , εxy )T и напряженийσ = (σxx , σyy , σxy )T состоит из трех групп уравнений: кинематических соотношений, определяющих уравнений и уравнений равновесия в областитела.

При отсутствии начальных напряжений в теле эта система уравнений27для плоского случая может быть записана в следующем матричном виде  εxx∂/∂x0u εyy  =  0∂/∂y  x ,uy2εxy∂/∂y ∂/∂x σxxd11σyy  = d21σxyd31∂/∂x00∂/∂yd13εxxd23   εyy  ,d332εxy σ∂/∂y  xx fσyy = x .∂/∂xfyσxyd12d22d32Записанная система трех матричных уравнений содержит вектор объемных сил с компонентами fx , fy , входящий в уравнение равновесия, матрицуупругих модулей с компонентами dij , связывающую напряжения и деформации в точке тела, а также две символические матрицы, состоящие изчастных производных по пространственным координатам.В матричном виде эта система уравнений имеет видε = Bu,σ = De,B T σ = f,(15)где D = (dij ) — симметричная матрица упругих модулей, B — матрица,состоящая из частных производных, f — вектор объемных сил.В задачах теории упругости различают плоско-напряженное и плоскодеформированное состояния.

С вычислительной точки зрения эти состояния будут отличаться только матрицей упругих модулей1 ν0E ,ν 10D=1 − ν20 0 (1 − ν)/21ν/(1 − ν)0E(1 − ν)ν/(1 − ν),10D=(1 + ν)(1 − 2ν)00(1 − 2ν)/(2(1 − ν))где E — модуль Юнга материала, ν — коэффициент Пуассона.Различают следующие граничные условия:1.

Кинематические граничные условия (условия Дирихле) u|Γ1 = ū. Например, часть поверхности закреплена.28(−1, 1) s(1, 1)r(1, −1)(−1, −1)Рис. 6: Четырехугольный билинейный элемент в локальной системе координат2. Силовые граничные условия (условия Неймана), когда на часть области действует сила. Действующие силы могут быть объемнымиf = (fx , fy ), поверхностными p = (px , py ) и сосредоточенными (узловыми) G = (gx , gy ).Система (15) эквивалентна задаче минимизации полной потенциальнойэнергии тела и сводится к задаче минимизации функционала энергии.Элементная матрица жесткости и элементный вектор правых частейдля плоской задачи с единичной толщиной элемента имеют вид ZZZfxpxeTee Te Tk =BN DBN dS, f =(N )dS +(N )d`, (16)fpyySeSe`eгде BN = BN — матрица производных от функций форм, S e — площадьэлемента, ` — ребро, на котором задана нагрузка.Для четырехугольного билинейного элемента (рис.

6) матрица BN зависит от координат и ее нельзя вынести за знак интеграла, как в случаелинейного треугольного элемента. Матрица системы имеет видXXKu = f, K =ke, f =f e + G,eeгде G — вектор узловых сил.Запишем формулу (16) для билинейных прямоугольных элементов. Длявычисления интегралов в (16) используем квадратурные формулы ГауссаZ1−1f (x)dx 'sXWi f (ti ).i=1Узлы t1 , t2 , . . . , ts и коэффициенты W1 , W2 , . . . , Ws подобраны так, чтоформула интегрирования точна для всех многочленов степени 2s−1. Тогда29элементная матрица (16) для билинейного элемента с использованием двухточек интегрирования для каждой координаты записывается в видеZTke =BNDBN dS =Se=2 X2XTWi Wj BN(ri , sj )DBN (ri , sj ) det |J(ri , sj )|,i=1 j=1∂∂xBN = B · N T =  0∂N1∂x= 0∂N1∂y0∂∂y0∂N1∂y∂N1∂x∂ ∂y∂∂xN10∂N2∂x0∂N2∂y0N10∂N2∂y∂N2∂xN20∂N3∂x0∂N3∂y0N2N300∂N4∂x∂N3∂y∂N3∂x0∂N4∂y0N30N400N4=∂N4 ∂y  .∂N4∂xПроизводные функций форм получают с помощью матриц Якоби∂Ns∂Ns ∂x −1  ∂r  ∂Ns  = J  ∂N  ,s∂y∂sPгде Pматрица Якоби вычисляется исходя из соотношений x =k Nk xk ,y = k Nk yk P ∂NkP ∂Nk∂x ∂yxk(ri , sj )yk(ri , sj )∂r∂r ∂r ∂r  kkJ(ri , sj ) = = .

(17)PP∂x ∂y∂Ni∂Nkxk(ri , sj )yi(ri , sj ).∂s ∂s∂s∂sikRРассмотрим вычисление интегралов вида ` N T · p d`, которые присутствуют в правой части уравнения (16). Пусть в качестве стороны интегрирования выступает сторона 2-3 прямоугольника. Эта сторона соответствуетзначению локальной координаты r = 1. Соответственно, функции формыбудут иметь видN1 = (1 − r)(1 − s)/4 = 0,N2 = (1 + r)(1 − s)/4 = (1 − s)/2,N3 = (1 + r)(1 + s)/4 = (1 + s)/2,30N4 = (1 − r)(1 + s)/4 = 0,Z 1Z` 1 TN · p ds =N T · p d` =N T | det(J)|ds =2 −1`−1 0000 000 0 px 0   Z 1 − sZ 1 px (1 − s) py ` 1``01 − s p xpy (1 − s) .

(18)=ds =ds = 0  py2 −1 1 + s2 −1 px (1 + s)2px  01 + spy (1 + s) py  00000000ZОсновные шаги программы. Используемые данные: n — число узлов; m — число элементов; массив mesh[n][2] = (xi , yj ) — сетка (координаты узлов) i, j = 1, n; массив elements[m][4] = (node1, node2, node3, node4) —список элементов (список узлов из которых состоит элемент); элементнаяматрица kij размера 8 × 8; элементный вектор правых частей b = (bi ),i = 1, 8; вектор правых частей f = (fi ), i = 1, 2n; orderGauss — порядок квадратуры; quadrGauss — массив размера orderGauss, содержащийточки интегрирования; weights — массив размера orderGauss, содержащий веса для точек интегрирования; u = (ui ), i = 1, 2n — вектор решения; us = t — граничное условие Дирихле; свойства материала — массивmaterial[m, 2] (первый индекс — номер элемента, второй индекс — материал: модуль Юнга E и коэффициент Пуассона ν; например, material[1, 1]содержит значение E первого элемента, material[20, 2] — содержит значение ν двадцатого элемента).Алгоритм1.

Считываем сетку, элементы, материал из файла — формируем массивы mesh, material, elements.2. В цикле по элементам: для i = 1 до m2.1 Определяем номера узлов текущего элемента: s1 , s2 , s3 , s4 . Длякаждого узла число переменных равно двум. Введем массив сномерами переменных numbers[8]Для j = 1 до 4 {sj = elements(i, j);numbers[2j − 1] = 2sj − 1;numbers[2j] = 2sj .

}312.2 Формируем элементную матрицу k и элементный вектор правыхчастей b (процедуру формирования см. ниже).2.3 Проводим процесс сборки: формируем глобальную матрицу жесткостиДля j = 1 до 8Для l = 1 до 8a[numbers[j], numbers[l]] = a[numbers[j], numbers[l]] + kj,l .2.4 Формируем глобальный вектор правых частейДля j = 1 до 8 fnumbers[j] = fnumbers[j] + bj .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. Выводим результат.6. Определение расчетных величин в элементах..Процедура формирования локальной матрицыи вектораR жесткостиTTDBNdV .правой части для задачи упругости (п. 2.2) k = V BN1. Определяем координаты узлов текущего элементаДля i = 1 до 4 {xi = mesh[si , 1]; yi = mesh[si , 2].}2. Определяем матрицу материала D текущего элемента – в зависимости от того, какое состояние плоско деформированное или плосконапряженное.3. В цикле по всем точкам интегрированияДля i = 1 до orderGaussДля j = 1 до orderGauss3.1 r0 = quadrGauss[i]; s0 = quadrGauss[j];3.2 Вычисляем матрицу Якоби J в точке (r0, s0) по формуле (17)3.3 Вычисляем определитель матрицы Якоби det J.323.4 Вычисляем обратную матрицу Якоби invJ.3.5 Вводим два массива для производных функций форм в ЛСК∂NiidN r[4] = ∂N∂r , dN s[4] = ∂s и вычисляем производные в точкеr0, s0dN r[1] = (s − 1)/4; dN s[1] = (r − 1)/4;dN r[2] = (1 − s)/4; dN s[2] = (−r − 1)/4;dN r[3] = (s + 1)/4; dN s[3] = (r + 1)/4;dN r[4] = (−s − 1)/4; dN s[4] = (1 − r)/4.3.6 Вводим два массива для производных функций форм в глобаль∂Niiной СК dN x[4] = ∂N∂x , dN y[4] = ∂y и вычисляем производныев точке r0, s0Для k=1 до 4 {dN x[k] = InvJ(1, 1)(dN r)[k] + InvJ(1, 2)(dN s)[k];dN y[k] = InvJ(2, 1)(dN r)[k] + InvJ(2, 2)(dN s)[k].}3.7 В точкеr0, s0 вычисляем матрицу BN = B · N = ∂N1∂N2∂N3∂N40000∂x∂x∂x ∂x∂N1∂N2∂N3∂N4 000= 0∂y∂y∂y∂y  =∂N1∂N1∂N2∂N2∂N3∂N3∂N4∂x∂y∂x∂y∂x∂y ∂ydN x[1]0.

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

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