Шабров Н.Н. - Метод конечных элементов в расчётах деталей тепловых двигателей (1061803), страница 15
Текст из файла (страница 15)
Структура массива 0И1 определена соотношением (4.94). Подпрограмма ЧКЛЛС формирует в соответствии с (4.94) и на основе сформированного ранее массива 0К1 матрицу Якоби и размещает ее в массиве К,1С. Параметрам ХРЕ, УРЕ и 7РЕ соответствуют одномерные массивы значений глобальных координат узловых точек Х!е!, 7!'! и Х!'! на основе соотношения (4.90).
В остальном работа программного модуля, реализующего вычисления по формуле (4.97), формально идентична рассмотренному ранее случаю для двухмерного плоского изопараметрического конечного элемента. Вычисление блока вектора узловых сил, эквивалентных действию распределенного в пределах шестигранного конечного элемента температурного поля, сводится к определению интеграла 1 1 1 Ге!~0!а — - — ) ) ~ В~'!Нео де1 ~ .1 ! с1~ !1!1 сЦ, (4.101) — 1 — 1 — 1 где матрица-столбец начальных деформаций е, дается соотношением (2.98). Матрица-столбец е, содержит в качестве скалярного множителя температуру, для которой можно по аналогии с (4.49) ввести интерполяционное соотношение вида (4.
102) Выполняя матричные преобразования в (4.101) и применяя одномерный метод численного интегрирования Гаусса, получим и л л Р,'~ = ~ ~ ~ с !ЗХ-1-2р!ОТ" де! 3 ~~, . ~ О~У,Н~, 1=! /=1 А=! !!а (4. 103) где Т!'! определяется выражением (4.102). 95 Фрагмент программного модуля допускает с учетом указанных изменений вычисление вектора узловых тепловых сил шестигранного элемента на основе представления (4.103) и размещение его в одномерном массиве ГЕ. Параметру ХРГ теперь следует присвоить значение, равное 3, т. е. числу степеней свободы в узле шестигранного конечного элемента. Рассмотрим вычисление вектора узловых сил, статически эквивалентных действию распределенной по поверхности шестигранного элемента нагрузки.
Для блока вектора Г~Я имеем Р. ~рй =- Е2Фа Рьь ь)а (К = 1, 2... т). (4.104) з(е) Рх Поверхностями шестигранного элемента являются грани этого элемента, для задания которых можно отказаться от использования трехмерных функций форм в (4.104), а воспользоваться независимой интерполяцией при помощи двухмерных функций форм. Итак, по аналогии с одномерным случаем (4.52) и (4.53) можно написать (4.105) г = Х У,Й, Ч)зт Где 7 — номер узловой точки на грани элемента; Ж вЂ” двухмерные функции формы, определенные в локальной системе координат независимо на каждой грани конечного элемента. 7 Поскольку в качестве распределенной нагрузки очень часто выступает нори мальное давление, то для — определения составляющих давления на глобальные оси г координат требуется вычислять направляющие кое, синуса внешней нормали к поверхности конечного зле~к мента.
Для того чтобы алгоРис. 4.5. К определениьо вектора едььнии- ритм вычисления направляной нормали к криволинейной поверх- ющих косинусов внешней ности трехмерного коььечьього элемента нормали не зависел от Грани элемента, следует установить правило обхода узловых точек, принадлежащих грани. Так, например, если смотреть на грань элемента со стороны внешней нормали, то обход узловых точек должен вестись ян прогив часовой стрелки так, как это сделано для плоского четырехугольного конечного элемента (см. рис.
4.3). В целях удобства обозрения рассмотрим грань шестигранника, для которой т1 = 1, и определим направляющие косинусы внешней нормали к ней (рнс. 4.5). Введем радиус-вектор точек, принадлежащих рассматриваемой грани, К =ХЕ,-1-УЕ, + гЕ2 (4. 106) где х, у, г даются соотношениями (4.105).
Выражение для орта внешней нормали будет иметь вид й, 2 Х К, и и=в 1й,-„,хи,,~ ' (4.107) где К,. и Й,„ — производные радиус-вектора К по локальным коордйнатам, Далее можно написать Ет Е2 Еб К,~х К,ч= х 2 у 2 г,~ =Ае1+Ве2+Сев, (4108) Хттг$гП ~й,ьх 1~,ч =ф 1~,В~ ~Р,~~ 1К,~'К,ч~ =~/ЕС вЂ” Р', (4.109) где А = у,;г, „ — у, „г, ; В = х,„,г, ~ — х, ~г,„,; (4.110) Е =- х ~+ у'в+ я'~,. 2 2 2 У, тт+2, Чт а =Х ВХ т1+У вВУ Ч+2 ~З Ч, (4.1 1 1) Тогда для направляющих косинусов внешней нормали можно написать: (4.112) Для элемента площади поверхности, заданной в форме (4.105) справедливо соотношение г)а= )Ц,, Х 1~, )с1~с1т1. (4.113) 97 4 шабров н.
н. соз(п, е,) = соз(п, е,) =- соз(п, е-,) =- А пе,=- 1й, хи,„!'! В С п'еЗ з Применяя к (4.104) одномерный метод численного интегрирования Гаусса и учитывая соотношения (4.112), можно последовательно написать ! ! л (е) Ррр = ~52Лу!!! — ! — ! п л Г!Р'! = РЮ2Л, А С А В С !1~ !1!1; (4.
114) Н;Н;. $!, ч! (4. 115) Для задания функции давления р в (4.115) можно воспользоваться ннтерполяционным соотношением на поверхности р=ХУ,В ч),, (4. 116) где рт — значение давления в узловой точке с номером у на соот. ветствующей поверхности. Глава Гз ИЗОПАРАМЕТРИЧЕСКИЕ КОНЕЧНЫЕ ЭЛЕМЕНТЫ ВЫСОКОГО ПОРЯДКА В ТЕОРИИ ТЕПЛОПРОВОДНОСТИ бл. ПЛОСКАЯ ДВУХМЕРНАЯ ЗАДАЧА Рассмотрим семейство изопараметрических четырехугольных конечных элементов первого и второго порядка применительно к решению плоской задачи стационарной теплопроводности.
Функции формы таких элементов и интерполяционпые соотношения для связи систем глобальных и локальных координат были установлены ранее при анализе плоской задачи теории упругости. Интерполяциопное соотношение для температуры в элементе будет иметь вид т!'! д, ч) = Е ж.д ~) т., (5.1) где Т вЂ” значение температуры в узле с номером а; Մ— функции формы, определенные в локальной системе координат соотношениями (4.28) или (4.29) (4.31); а = 1, 2, ..., т. Выражение для блока матрицы градиентов запишется так: д д.
Л'- (5.2) В," ,= д Д ~а Здесь преобразование производных функции форм по глобальным координатам в производные по локальным координатам осущест- вляется в соответствии с (4.38) — (4.40). эв Введем обозначение д — Л' а (5.3) — У„ ду а где Ь, с„— элементы массива ВЕ, определенного в (4,45). Тогда на основе соотношения (3.28) для элемента матрицы теплопроводности в данном случае можно написать ! ! Р('„!а = ~ ) (Л,Ь,,Ьа+ Х„с„са) с(е1 Л сй с(!), (5.4) — 1 †! где )!,, и Х„ — коэффициенты теплопроводности по осям х и у. Применение к (5.4) одномерного метода численного интегрирования Гаусса дает л л Рса(! =,Уд ~'д (РихЬаЬа + Хуст(!) сне( Л $ц ч НсНу.
(5.5) с=!1=! Из (5.5) видно, что выражение для элемента матрицы теплопроводности четырехугольного изопараметрического конечного элемента, вычисленное в одной точке интегрирования Гаусса, совпадает с точностью до скалярного множителя с выражением (3.29). Следовательно, для программирования выражения (5.5) можно использовать фрагмент программного модуля, приведенный на стр. 62, который должен быть окаймлен двумя циклами интегрирования по Гауссу с добавлением стандартного для изопараметрических элементов фрагмента в виде подпрограмм РИРМ1, РКГ)М1, РААС, М1МЪ" и МАТВЕЕ.
Ниже приведен фрагмент программного модуля для вычисления в соответствии с (5.5) элементов нижней симметричной части матрицы теплопроводности и размещения ее в виде одномерного массива ЬЕ. САВЕ ТАВГ.Е(ЯСт) ЛО 10 И=1,(Ча сАеь АЙАГ)зз(и,х,н1) Г!О Ы .)а=1,(Ча сАГ.ь ы()зз(,) ад,н,)) 2 СА1.1 РИЕЫ1(ЯРЕ,Х,У,Х!,У1,РИ!) 3 СА1.1. Р1(Э(Ч1( (РЕ,Х,У,Х1,У1,Г>(Ч1) 4 СА?Л РМЗАС(Ь!РЕ,ХРЕ,УРЕ,В!)Ч1ЯЛС) СА1.1. М11ЧЧ (КЛС,КЭМ,1)ЕТ,1.1,1.2) САГ.Г. )ЧЫВЕЕ((ЧВМ,ЯРЕ,Г!)Ч1,КЛС) 1 ЭЕТ=ЭЕТэН1эНЛ 00 20 11.=1,ЯРЕ 4ф Ь.З.
ТРЕХМЕРНАЯ ЗАДАЧА Рассмотрим трехмерные шестигранные изопараметрические конечные элементы первого и второго порядка применительно к трехмерной задаче теории стационарной теплопроводности. Функции формы и необходимые интерполяционные соотношения для этих элементов были установлены при анализе соответствующей задачи теории упругости. Интерполяционное соотношение для температуры в элементе будет Т (а, 11, ь) = Е Ла($~ ть ь)Ти, а (5.2!) где Ти — значение температуры в узле с номером и; ӄ— трехмерные функции формы, определенные в локальной системе координат соотношениями (4.82) или (4.83) — (4.86). ! Для блока матрицы градиентов можно записать — Л"„ — Л!а д — Л! дг а Ва (5.22) Преобразование производных функций форм по глобальным координатам. в (5.22) в производные по локальным координатам осуществляется в соответствии с (4,92) — (4.94).
Введем обозначение — Л!„ д д Ва' (5.23) где Ь„, си, 0а — элементы массива ВВ, определенного в (4.99). Нааа основе (3.28) для элемента матрицы теплопроводности шестигранного конечного элемента с учетом (5.23) можно записать ! ! ! Рсар = ~ ~ ~ (~хна(!р + ~~усагр ! Маг!р) !(е( ~ !)~ !(г! Й~ (5.24) †! †! †! где матрица Якоби 3 дается соотношением (4.94). Численное интегрирование в (5.24) по методу Гаусса дает а и а Рсйр = Е Е Е (~х~?а!!р+ Хусаср+ Хддас(р) СЫ Л ~ . я, !,НГН!Н!, ° ! — ! !.=! ! — 1 (5.25) Сопоставляя выражения (5.25) и (5.5) с точки зрения программирования, видно, что значения этих выражений, вычисленные 103 для одной точки интегрирования, программируются одним и тем же оператором с меткой 30 в программном модуле, приведенном на стр. 99.
Параметру МВМ должно соответствовать значение, равное 3. Далее очевидно, что для структуры программного модуля, реализующего вычисчения по формуле (5.25), будет идентична структура программного модуля, реализующего вычисления по формуле (5.5).
Фактическое отличие сводится к замене операторов с метками 2, 3, 4 последовательностью операторов, приведенных на стр. 94. Для элемента матрицы поглощения в соответствии с (3.15) можно написать ! ! ! Рфв= ~ ~ ~ КУ~Увде130~!1!1д!.. †! †! ! После численного интегрирования в (5.26) методом Гаусса получим л л а Р~а'~,в — Е Е ~ КУиУвде1 .1 ~~ „! Н!Н~Н!,. (5.27) с=! !=! й=! Для задания функции поглощения й можно воспользоваться интерполяционным соотношением ~6 т1 ~)=~У 6») ~)Л, а (5.26) (5.28) где ߄— значение функции Я в узле с номером я. Матрица конвекции шестигранного конечного элемента определяется в соответствии с (3.14) интегрированием заданной на поверхности функции. В связи с этим здесь также разумно будет воспользоваться двухмерными функциями формы для интерполяции поверхности трехмерного конечного элемента при помощи соотношений (4.
105). Тогда для элемента матрицы копвекции можно написать ! ! РД =- ~ ~ЬУ,,У~ ~/ Еб — Е'с5д0, (5.29) †! †! где У, ӄ— двухмерные функции формы, определенные независимо на поверхности конечного элемента, через которую осуществляется конвекция. Функции Е, 6 и Е даются соотношениями (4. 111) . Применяя одномерный метод численного интегрирования Гаусса к (5.29), получим л и Р!т! = ~ „~; ЬУтУь) Ег! — Е'; „,Н<Н;.
(5.30) ! !/ ! Для зада«ия функции Ь на поверхности конвекции конечного элемспта можно воспользоваться интерполяционным соотношением типа (4.116). !оа Легко получить в соответствии с (3.16) — (3.18) следующие формулы численного интегрирования для компонент векторов узловых тепловых сил: (5.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
