Сагдеева Ю.А., Копысов С.П., Новиков А.К. - Ввеление в метод конечных элементов (1061801), страница 4
Текст из файла (страница 4)
Функция f (x, y) в левой части равенства представляет собой функцию, выраженную в глобальной системе координат,18тогда как f (x(r, s), y(r, s)) соответствует функции, представленной в локальной системе координат (x, y — глобальные (старые) координаты, r, s— локальные (новые) координаты). Матрица Якоби имеет вид ∂x ∂y ∂r∂rJ = ∂x∂y .∂s∂sДля треугольного элемента наиболее распространенной является бариоцентрическая система координат или L-координаты. Каждая координата представляет собой отношение расстояния от выбранной точки треугольника до одной из его сторон s к высоте h, опущенной на эту сторонуиз противолежащей вершины (рис. 3). Величины Li изменяются в пределахот нуля до единицы.L-координаты треугольника удовлетворяют соотношениюL1 + L2 + L3 = 1.Уравнения этого типа следовало ожидать, потому что три координаты вдвумерном случае не могут быть независимыми.
Местоположение произвольной точки может быть полностью описано с помощью только двухкоординат.L-координаты обладают следующими свойствами. Координатные переменные L1 , L2 и L3 представляют собой функции формы для треугольногоэлемента: Ni = L1 , Nj = L2 , Nk = L3 (Li = 1 в узле с номером i и равна 0в узлах j и k).Координаты произвольной точки C в декартовой системе выражаютсячерез L-координаты следующим образомx = L1 xi + L2 xj + L3 xk ,y = L 1 yi + L 2 yj + L 3 yk .Преимуществом использования L-координат является существованиеинтегральных формул, которые упрощают вычисление интегралов вдольсторон элемента и по его площади:Za!b!(L1 )a (L2 )b d` =`,(a + b + 1)!`Za!b!c!(L1 )a (L2 )b (L3 )c dS =2S,(a+b+ c + 2)!S∆19u2kkxujuiiyu2k−1jxu2ijiu2ju2j−1u2i−1Рис. 4: Степени свободы векторных величин в узлах элементовгде S — площадь элемента.
Использование этих соотношения может бытьпроиллюстрировано при вычислении интеграла видаZZ2S1!1!0!2S == S/12.Ni Nj dS =(L1 )1 (L2 )1 (L3 )0 dS =(1 + 1 + 0 + 2)!4!S∆S∆7Аппроксимация векторных величинДо сих пор рассматривалась интерполяция скалярных величин, такихкак, например, температура или давление. Скалярная величина представлялась как линейная комбинация функций формXue =u i Ni .Векторные величины (например, перемещение) принято разделять на компоненты или степени свободы. В зависимости от размерности задач каждый узел содержит одну, две или три неизвестных (рис.
4).В одномерном случае представление векторов и скалярных величинсовпадают ui.ue = Ni ui + Nj uj = (Ni Nj )ujВ двумерном случае отдельно аппроксимируются горизонтальная ивертикальная компоненты с помощью одних и тех же базисных функций.Перемещения вдоль осей x, y обозначим через ux , uy . Для треугольного20элементаux = Ni u2i−1 + Nj u2j−1 + Nk u2k−1 ,uy = Ni u2i + Nj u2j + Nk u2k .или в матричном виде: uxNi=uy080NiNj00NjNk0u2i−1 u2i 0 u2j−1 .Nk u2j u2k−1 u2kПрименение треугольных конечныхэлементов в задаче теплопроводностиПостановка задачи. Стационарное уравнение теплопроводности имеетвид−∂∂T∂∂T∂∂T(Dxx)−(Dyy)−(Dzz) + Q = 0,∂x∂x∂y∂y∂z∂z(7)где T — температура, Dxx , Dyy , Dzz – коэффициенты теплопроводности(Вт/(м град)), Q — плотность тепловых источников внутри тела (Вт/м3 ),эта величина считается положительной, если тепло подводится к телу.Граничные условия могут быть трех типов:1.
На границе тела (или ее части) задана температура, которая можетбыть функцией координатΓ1 : T = T̄ (x).(8)2. На границе задан конвективный теплообмен. Конвективный теплообмен — процесс переноса тепла, происходящий в движущихся текучихсредах (жидкостях либо газах) (теплообмен с окружающей средой).Конвективный теплообмен характеризуется величиной h(T −T∞ ), h —коэффициент теплообмена Вт/(м2 · град), T∞ — температура окружающей среды.3. На границе задан тепловой поток q (то есть извне подводится или отводится тепло).
Тепловой поток q положителен, если тепло отводитсяот тела.21Граничные условия п. 2 и 3 записываются с помощью смешанного условияΓ2 : Dxx∂T∂T∂Tcos α + Dyycos β + Dzzcos γ + h(T − T∞ ) + q = 0,∂x∂y∂z(9)где (cos α, cos β, cos γ) — единичный вектор внешней нормали к поверхности. Если конвективный обмен отсутствует, то есть отсутствует обмен тепла с окружающей средой и поток тепла равен нулю, условие на границеприобретает вид ∂T /∂n = 0 — условие теплоизолированности границы.Решение уравнения (7) с граничными условиями (8)–(9) эквивалентноотысканию минимума функционалаZZZZF (T ) =(∇T )T D∇T dV + 2T qdΓ +h(T − T∞ )2 dΓ + 2QT dV.Γ2VΓ2VРазобьем область V на конечные элементы Ve , e = 1, m, и аппроксимиnPруем решение T конечной суммой T =Ni Ti , где n — число узлов, Ni —i=1функции формы.
Записывая условие минимума функционалалучим СЛАУ относительно узловых значений:XXkeu =feeK=Xke ,F =eXf e.eЭлементная матрица системы имеет видZZke =(B e )T De B e dV +h(N e )T N e dΓ = k1 + k2 ,Vefe = −ZΓ2= 0, по-eилиKu = F,∂F∂Ti(10)Γ2q(N e )T dΓ +ZQ(N e )T dV +Zh(N e )T T∞ dΓ.(11)Γ2VРассмотрим треугольный элемент с тремя узлами i, j, k.
Функции фор1(as + bs x + cs y), s = i, j, k,мы линейного треугольного элемента Ns = 2A(где as , bs , cs определяются соотношением (6)). Температура на элементеопределяется формулой TiT e = Ni Nj Nk T j .Tk22Запишем матрицу градиентов от функций формы1 bi bj bkBe =,2A ci cj ckматрица свойствD=Dxx00.DyyТеперь можно вычислить матрицу теплопроводности элемента (10) (элементную матрицу жесткости). Первое слагаемое k1 принимает видZZb i ci 1Dxx0bi bj bke T e eb j cjdV.(B ) D B dV =0Dyyci cj ck4A2 V eVeb k ckПредполагая толщину элемента единичной, заменим dV на dA. Верхнийиндекс e опускаем. Подынтегральное выражение постоянно и может бытьвынесено за знак интеграла:ZZTTB DBdV = B DBdA = AB T DB.VeAВычисляя произведение матриц, имеемb b bi bj bi bkci ci ci cj ci ckDxx i iDyy bj bi bj bj bj bk +cj ci cj cj cj ck .k1 =(12)4A4Abk bi bk bj bk bkck ci ck cj ck ckRВторой интеграл S hN T N dS должен быть вычислен по поверхностиZZNi Ni Ni Nj Ni Nkk2 =hN T N dS = h Nj Ni Nj Nj Nj Nk dS.(13)SSNk Ni Nk Nj Nk NkФункции формы зависят от x и y, поэтому произведения вида Ni Nj немогут быть вынесены за знак интеграла.
Кроме того, значение интегралазависит от того, на какой поверхности наблюдается конвективный теплообмен. Если, например, конвекции подвержена сторона между узлами i и j,то Nk равно нулю вдоль этой стороны и интеграл сводится к следующемувыражениюZZNi Ni Ni Nj 0(14)hN T N dS = h Nj Ni Nj Nj 0 dS.SS000234516342321Рис. 5: Пример сетки из треугольных элементов. Перечислены номера узлов и номера элементовномерэлемента1234узелi1122узелj5263узелk4556Таблица 1: Списки узлов из которых состоит элемент (массив elements)Если любая из двух других сторон подвержена конвекции, то расположение отличных от нуля членов в (13) будет иным, чем в (14).Вычисление произведений от функций форм не представит труда, еслиприменить L-координаты и интегральные формулы для них.
Три интеграла в выражении для вектора нагрузки элемента также легко вычисляются,если воспользоваться L-координатами.Программная реализация. Используемые данные: n — число узлов,m — число элементов; сетка (координаты узлов) mesh[n][2]: mesh[i][1] — координата xi , mesh[i][2] — координата yi , i, j = 1, n; массив elements[m][3] =(node1, node2, node3) — список элементов (списки узлов из которых состоитэлемент), для сетки на рисунке 5 массив elements имеет вид: elements[1, 1] =1, elements[1, 2] = 5, elements[1, 3] = 4 (см. Табл.
8); элементная матрицаkij размера 3 × 3; элементный вектор правых частей b = (bi ), i = 1, 3;вектор правых частей f = (fi ), i = 1, n; u = (ui ), i = 1, n — вектор решения, us = t — граничное условие Дирихле, свойства материала — массивmaterial[m, 2] (первый индекс — номер элемента, второй индекс коэффициент теплопроводности Dxx , Dyy ), так material[1, 2] есть значение Dyyпервого элемента, material[20, 1] есть значение Dxx двадцатого элемента.Толщина элемента предполагается единичной.24Алгоритм1. Считываем сетку, элементы, материал из файла — формируем массивы mesh, material, elements.2. В цикле по элементам: для i = 1 до m2.1 Определяем номера узлов текущего элемента:Для j = 1 до 3 {sj = elements(i, j)}.s1 , s2 , s3 — номера узлов текущего элемента.2.2 Формируем элементную матрицу k = (kij ): k11 , .
. . , k33 (формирование матрицы для задачи теплопроводности см. ниже).2.3 Проводим процесс сборки: формируем глобальную матрицу жесткости.Для j = 1 до 3Для r = 1 до 3 {asj ,sr = asj ,sr + kj,r }.2.4 Формируем элементный вектор правых частей b = bi .2.5 Проводим процесс сборки — формируем глобальный вектор правых частейДля j = 1 до 3 {fsj = fsj + 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. Решаем систему Ku = f .5. Выводим результат u.6.
Определение расчетных величин в элементах.Процедура формирования локальной матрицы жесткости для задачитеплопроводности (п.2.2 алгоритма, формулы (10), (11), (12), (14)).1. Определяем координаты узлов текущего элементаДля i = 1 до 3 {xi = mesh[si , 1]; yi = mesh[si , 2]}.252.
Формируем 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. Формируем матрицу k1 по формуле (12).5.
Если на одной из сторон задан конвективный теплообмен, напримерна стороне 12, то вычисляем длину стороны L12 и матрицуZ2 1 0hL12 1 2 0 .k2 =hN T N dS =6S120 0 06. Формируем матрицу k = k1 + k2 .7. Если на одной из границ задан тепловой поток, его надо учесть вправой части. Пусть напримерпоток задан на стороне 23. Вычисляpем длину стороны L23 = (x3 − x2 )2 + (y3 − y2 )2 . К вектору правыхчастей добавляем вектор Z0N1N2 dS = qL23 1 .q2S1N38. Если на одной из границ задан конвективный теплообмен, его надоучесть в правой части.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
