Шабров Н.Н. - Метод конечных элементов в расчётах деталей тепловых двигателей (1061803), страница 14
Текст из файла (страница 14)
Сд $! ч; (4.75) Вычисление значений радиуса и температуры в точке интегрирования достигается подпрограммой-функцией ЧА1.УЕ, работа которой обьяснялась выше. Формальному параметру ХРЕ теперь соответствует одномерный массив значений радиусов узловых точек элемента. Приведенный фрагмент программного модуля включает вычисление вектора узловых тепловых сил конечного элемента Г,!! (е) и . размещение его в одномерном массиве ГЕ .
Б этом фрагменте следует обратить внимание на то, что вычисление вектора узловых тепловых сил в соответствии с выражением (4.75) и вычисление матрицы жесткости элемента на основе (4.70) и (4.71) объединены общими циклами интегрирования по методу Гаусса. Это неявно предполагает, что порядки интегрирования в (4.70), (4.71) и (4.75) совпадают. Очевидно, что такое совпадение порядков не является обязательным. В связи с этим для вычисления вектора тепловых сил и матрицы жесткости элемента требуется организовать раздельные циклы интегрирования, что приводит к существенному усложнению логической структуры программного модуля. Выражение для блока вектора узловых сил, статически эквивалентных действию распределенной по поверхности тороидального элемента нагрузки, будет иметь вид Гр',„, — — 2й Е2Л~,у ~ 1/ й (сс — — 1т 2, ..., Ш).
(4.76) !. Рг1 87 г = ~~ У (~) г ; (4.77) (4.78) где у — номер узловой точки на стороне элемента; г,„и г глобальные координаты узловых точек стороны элемента. Тогда по аналогии с (4.60) для вектора узловых снл элемента (4.76) можно написать следующее выражение: (4.79) где для р справедливо интерполяционное соотношение (4.62). Применяя к (4.79) одномерный метод численного интегрирования Гаусса, получим (4.80) Здесь $, — точки интегрирования, а Н, — весовые коэффициенты Гаусса. Обсуждение вопросов программирования выражения (4.80) будет сделано позже. В табл. 4.1, 4.2 приведены результаты расчета круглой плиты постоянной толщины под действием температурного поля, изменяющегося по квадратичному закону вдоль радиуса.
При одинаковом числе элементов результат, полученный с использованием конечных элементов второго порядка, дает хорошие совпадения с точным решением даже для напряжений. Примечательно, что практически точно удовлетворяются краевые условия на внешнем контуре плиты.
4.4. ТРЕХМЕРНАЯ ЗАДАЧА Рассмотрим семейство трехмерных изопараметрических шестигранных конечных элементов первого и второго порядка, для которых интерполяционные полиномы являются соответственно линейными и квадратичными функциями локальных координат $, 90 Будем считать, что распределенная по поверхности нагрузка представляет собой нормальное давление р. Следовательно, все рассуждения относительно определения направляющих косинусов внешней нормали к контуру элемента, проведенные в и. 4.2, остаются справедливыми и здесь. Для задания контура элемента, на котором задана распределенная нагрузка, воспользуемся одномерными функциями формы и в соответствии с (4.52) и (4.53) напишем т) и ь (рис.
4.4). Пределы изменения локальных координат для всех элементов составляют — 1 -$, т), ~ ~ 1. (4.81) Шестигранный элемент первого порядка имеет восемь узловых точек, которые расположены в вершинах элементов (рис. 4.4, а). Функции формы такого элемента даются соотношениями Ж;6, Ч, 3= — (1+Бд(1+Чт)д(1+Кд, (4.82) где й;, т);, ~; — локальные координаты узловой точки с номером 1, ( = 1, 2, „ 8. г Рис.
4,4. Линейный (а) и квадратичный (б) иэопараметрические трехмерные конечные элементы в глобальной и локальной системах координат Шестигранный элемент второго порядка имеет двадцать узловых точек, первые восемь из которых расположены в вершинах элемента, а остальные — на средине ребер (рис. 4.4, б). Функции формы такого элемента даются соотношениями: Лг~б т), ') = — (1+Бд(1+Чт)д(1+1~д(Б;+ т)Ч;+Б; — 2) (Е= — 1,2, ...,8); (4.83) Лг~й т), .".) = 4 (1 — т)')(1+Му)(1+Гу) (1'=9, 11, 17, 19); 6 т) ) == — (1 — )(1+ т)тЬ)(1+ Е~ ) Т (е, ), 1) == — ',1 — ."') (1 + ".еь,) (1 -' Ид (4.84) (й = 10, 12, 18, 20); (4.85) (1 = 13, 14, 15, 16). (4.86) 91 Интерполяционные соотношения для перемещений в элементе будут иметь вид (4.87) где пг — число узлов элемента; и„, п„, ы„ †компонен вектора перемещений в узлах.
Аналогичным образом строится интерполяционное соотношение для формы элемента, которое можно написать т х (е> а е,у„е а=( хе (4.88) где х, д, г„— глобальные координаты узловой точки с номером (х. Соответствие между глобальными и локальными координатами точек элемента можно представить и в другой форме: хД, т(, Д=М(~Д, т1, 1)Х(~; (4.89) Здесь М(е1 — матрица функций форм элемента, а Х(е), т'(е>, Х(е> векторы узловых значений глобальных координат, т.
е. х, Х2 Х(о = у(е) (е) ъ Х (4.90) Выражение для блока матрицы градиентов имеет вид Р~й = РЕзЛ(а~ (4.91) где Р— матричный дифференциальный оператор — дается выражением (2.90); ܄— функции формы элемента — даются соотношениями (4.83) или (4.84) — (4,8б); а = 1, 2, ..., т. Компонентами блока матрицы градиентов (4.91) являются производные по глобальным координатам от функций форм, заданных в локальной системе координат, Следовательно, здесь также требуется выполнить весь комплекс преобразований, устанавливающих связь между производными функций форм по гло- 92 В соответствии с (4 44) введем обозначения д А~а д (4.98) Тогда произведение В„"НВ~~' в (4.97) дает результат, который с точностью до множителя 1/36 (1'('1)' совпадает с результатом в аналогичном произведении в (2.94) для трехмерного симплекс- элемента.
Это обстоятельство дает возможность воспользоваться без изменения подпрограммой ЯТ1ГГ для вычисления блока ма- трицы жесткости шестигранного изопараметрического конечного элемента в одной точке интегрирования Гаусса. Двухмерный массив ВЕ в данном случае также предназначен для хранения ненулевых элементов матрицы градиентов В(') и определяется следующим образом Ь1...(а ... д„, ВЕ = с1 ... с„... с,„ (4.99) ~11 ° ° ° ~(а ° ° 11т Параметру РЕТ присваивается значение 1)ЕТ = г(е1 Л (~п т)~, ~1) Н;Н Нд. (4.100) Здесь также используется схема хранения нижней симметрич- ной части матрицы жесткости К(') в виде одномерного массива, при которой можно без изменения воспользоваться подпрограммой сортировки ГКЧ5Е элементов каждого блока К$ в одномерный массив 5Е.
Очевидно, что структура программного модуля, реализующего вычисления по формуле (4.97), будет идентична структуре при- веденного фрагмента программного модуля и реализующего вы- числения по формуле (4.43). Фактическое отличие заключается в добавлении нового цикла интегрирования по локальной коорди- нате ~ и замене подпрограмм РКГМ1, РЕВЫ и РААС подпро- граммами ЧМГИ1, ЪК1ЭХ1 и УМЛАС (см.
приложение). Таким образом, все исправления в упомянутом программном модуле сводятся к замене операторов с метками 1, 2 и 3 нижеследующими операторами: 00 10 КО=1,ИО СА!Л. ЯЛ()Ю (КС,Е,НК) САП. У11РХ! (КРЕ,Х,У,Х,Х1Л1,71,РИ1) СА1 1- УМРИ! (МРЕ,Х,У,2,Х1,71,2!,ЛИ1) СА1Л УМ)АС (ЯРЕ,ХРЕ,"1'РЕ,ЕРЕ,1ЭМ!,Й)С) Подпрограмма ЧМГК1 формирует для каждой выбранной точки интегрирования при значениях Х = $;, У = Ч1г — Я, = Ь1, 94 одномерный массив РИ1, компонентами которого являются значения функций форм конечного элемента.
При ЯРЕ = 8 формируется массив ЕЫ для конечного элемента первого порядка в соответствии с соотношениями (4.82). При ЯРŠ— -- 20 формируется массив ГЫ для конечного элемента второго порядка в соответствии с соотношениями (4.83) — (4.8б). Параметру ЯРЕ, как и прежде, отвечает число узлов конечного элемента. Формальные параметры Х1, У1, Х1 представляют собой одномерные массивы локальных координат узловых точек соответственно по $, Ч, ь. Подпрограмма ЧХ1ЭМ1 формирует для точки (51, 1~;, ь!) интегрирования двухмерный массив 1ЭЫ1, столбцами которого в зависимости от значения параметра ХРЕ являются производные по локальным координатам функций форм шестигранного конечного элемента первого или второго порядка.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
