Сегерлинд Л. Дж. - Применение метода конечных элементов (1051193), страница 20
Текст из файла (страница 20)
Пустые места в [К~ означают ну. левые коэффициенты. Значение Т, известно (150'С), так что система уравнений должна быть модифицирована перед решением. Эта модификация преобразует столбец правых частей к виду (Р)т 18700 7650 1200 1200 1200 1000) После решения системы имеем (Т)т 1150 826 59 486 442 4261 Теоретические значения температуры 121 следующие: (Ттсо „] т =1150 89>9 62,8 50,6 45,2 43>31. Результаты, полученные-по методу конечных элементов, доста.
точно хорошо согласуются с истинными значениями, если учесть, что было проведено разбиение области на одинаковые элементы, 142 Глана а Решение по методу конечных элементов можно было бы улучшить, если использовать более короткие элементы вблизи стены, в кото. рую заделан стержень.
В предыдущем примере площадь поперечного сечения стержня была постоянной. Этого требовать необязательно. Площадь может быть также функцией длины элемента. Если площадь элемен- Фнг. 8 2. Одномерный элемент с переменной площадью сечения. та меняется по длине, то,матрица элемента должна быть преобразована. Обсудим теперь это преобразование. Рассмотрим элемент, изображенный на фиг. 8.2. Его площадь поперечного сечения изменяется от А» на левом конце до А; на правом. Если ограничиться случаем линейного изменения площади, можно сразу записать соотношение для площади А =»т'»А»+ У»А;, где У» и У» — линейные функции формы. Следующий шаг — определение матриц элемента.
Используя для вычисления 1й»о) первую часть формулы (8.15), получаем (В)г (0) (В)»й»= — ~ ' ~ ~ А»(х. (8.28) У о Так как величина А в этой формуле не постоянна вдоль х, ее нельзя вынести за знак интеграла. Подстановка выражения (8.27) в (8.28) дает — 1 — — А, + — А» дх. о Выполняя интегрирование, получаем — — А,+ — А,, Леренос тепла еа счет теплопроводности и конвекции или де 1 1 А+А! (8,29) Так как 1Ае+Ат)/2 равно средней площади элемента, выражение (8,29) может быть записано как (8.30) Р=М,Р,+1т' Ртэ (8.31) с учетом которого имеем (п1)т (п1),1с ь ( ~ ~ т т Р (х ,) 1Ф$у! Ф1 3 о (8.3лу Первый коэффициент в (8.32) после подстановки У» и У; и вычисления интеграла будет равен ~И,'Рс(х =~ (ЩР, + М;.И;РД Ых= — — (3Р, + Рт).
(8.33) о о Поверхностные интегралы представляются соотношениями й (1У)тРУ),15 — —" ( '+ «) ( + ~) (8.34) " ~(Р~+Ь ж+3Р 1 ЬТ (У1~ сЖ= (8.35) Сложив матрицы (8.30) и (8.34), можно получить (94]. Соот. ношение (8.35) дает часть Д<е>). Из формул (8.34) и (8.35) ясно, что их нельзя получить из (8.17) или (8.21) простой подстановкой среднего периметра в случае элемента конической формы. где А — средняя площадь.
Формула (8.30) совпадает с (8.15) о точностью до замены площади,на ее среднее значение. Интегралы по боковой поверхности могут быть выражены аналогичным образом. Периметр можно записать в виде соотношения [46 Глаза 8 Если предположить, что конвективный теплообмен имеет место на поверхности стороны элемента между узлами 1 и 1, то в точках этой поверхности Мь=Е4=0 и соотношение (8.44) примет вид Е,Е, Е,1 О ь [Цг [А~[,(5 3 УО (8.45) где Ю=ЫЯ, причем толщина ! предполагается единичной.
Два типа произведений входят в формулу (8.45): квадрат величины Е,' или Ц' и перекрестное произведение Е|Еь Начнем с квадратных членов: 2 1 О й ~~ЬГ~г [А![~КЯ= — ' 1 2 О Мп О О О (8.46а) Аналогичные соотношения получаются для стороны между узла- ми)и яг й[И[г [А!) (Я= — "'",' ~м (8.46б) и для стороны между узлами Й и 1: (8.46в) Три интеграла в выражении для вектора нагрузки элемента также легко вычисляются, если воспользоваться Е-координатами.
Начнем с интеграла ) [йг!тЯЫУ; предположим, что величина Я = —" ° (2+0+ !)! '~ 3 ~и ~п гдето — длина стороны между узлами 1 и 1. Интегрирование перекрестного произведения дает [![Мп,йц ЕгЕз~~= ([+!+ц! а Мсс Интегралы ) ЕзШ и ) ЦЫ равны между собой. Подставляя Йп Ж~ полученные результаты в формулу (8,45), получаем Перенос текла еа счет тенлоироеодности и конеекции постоянна внутри элемента. Тогда будем иметь д ~рц)г,(у д ~ У (8.47) 0 (8.49а) ~~ (цт дс (8.496) (8.49в) Величина йт„~()У) (8 идентична (8.49а) — (8.49в) с учетом замены д на йТ . Если тепловой поток или конвективный теплообмен наблюда.
ются на двух сторонах элемента, то поверхностный интеграл заменяется суммой интегралов по каждой из сторон. То же самое относится к интегралу (8.43). В большинстве задач о переносе тепла интерес представляют значения температуры в узловых точках. Иногда бывает необходи- Таким образом, тепло, генерируемое в элементе, распределяется поровну по трем узлам.
Интегралы ~ [Лт1тдс(5 и ~ [)у1тйт„с(5 за Зт е писываются в одинаковой форме: Ж1 Ф1 с(5, (8.48) Лтн так что только один из них необходимо вычислить. Поскольку интеграл (8.48) поверхностный, его можно рассмотреть так же, как интеграл (8.43). Результаты зависят от того, на какой из сторон элемента происходит конвективный теплообмен, характеризуемый величиной Ь, или приток тепла за счет теплового потока д. Предполагая д постоянным по поверхности элемента, получим для интеграла (8.48) три следующие формы записи; 148 Гааза д дт дх — та (8.50) дт ду Численный пример, который приводится ниже, иллюстрирует применение полученных выше соотношений. Пример 60.
Ниже изображен элемент, который использован для дискретизации сплошной среды. По двум поверхностям этого элемента происходит конвективный теплообмен. Указаны размеры элемента и физические характеристики. Требуется составить матрицы элемента, предполагая его толщину единичной.
'с Фм'Кз К Кут дд Вт((ем -Ю К задаче 60. Рассмотрии матрицу теплоцровод~ности элемента: Д(е)! 4А ь,ь, ь,ь; ь,ь ь(ь, ь ь ь ь ь(ьа Ь|ьа ьаьа С;С; С(С( С(С1 С1С1 с(са сзса "[ мо определять градиенты температуры. После того как определе- ны узловые значения, градиенты температуры находятся с помо- щью соотношения Перенос тепла ва счет теплопроводности и конвекции Константы Ь и с вычисляются по формулам Ь,=1; — )ел=Π— 4= — 4, Ь, =1'л — )т; =4 — 3 =1, Ь =['; — 1' =3 — 0=3, с,=Մ— Хг — — 6 — 7= — 1, с =Х,— Х„=З вЂ” 6= — 3, с„=Х; — Х;=7 — 3=4, 1 Х; Х )т 1 3 3 1 7 0 1 6 4 2А= =13, А= 2 Длины сторон Я!и и Хы равны ПОдетаНОВКа ПОЛУЧЕННЫХ ЧИСЛОВЫХ ЗНаЧЕНИй В [Ь1е') даЕт 16 — 4 [л1е1) = — — 4 1 36 26 — 12 3 1 3 — 4 3 9 — 12 — 4 — 12 16 6 (4,12) 6 или 6,32 0 3,16 0 8,24 4,12 3,16 4,12 14,56 17 — 1 à — 1 10 — 16 — 9 [!С(е)[ 26 24,88 — 1,15 21,09 — 1,!5 18,41 — 6,96 .
21,09 — 6,96 40,98 Вектор нагрузки элемента ([<е~) представляет собой сумму двух поверхностных интегралов по каждой из сторон, где происходит 21и — — )/(Хг — Хи)е+()ст — )тл)е )У(7 — 6)е+(Π— 4)т=)/17=4,12 см„ е = У~к — х) ч ~У вЂ” У,с=УГе — ЗГГч4 — 3)'=УГ0=316 Перенос. тепла еа счет теплопроеодности и конеекции -О О О О)ссе!)т !)сс) Дс полы 12 О 1 2 1 О 1 1 2 (8.55) Рй)„)р О (8.56) О 1 1 1 1~ц ) = 3 (8.57) Для интеграла (8.55) существуют три другие формы записи, по одной на каждую из оставшихся сторон.
В каждой из них значения коэффициентов на главной диагонали равны двум и значения ненулевых коэффициентов вне главной диагонали равны единице. Коэффициенты в строках и столбцах, соответствующих узлам, расположенным вне рассматриваемой поверхности, равны нулю. Для интеграла (8.57) тоже существуют три другие формы записи.
Нулевой коэффициент находится в строке, соответствующей узлу вне рассматриваемой поверхности. Юоп — площадь поверхности, содержащая узлы с, 1, й и т. д. 8.5. Преобразования координат При выводе уравнения (8.1) используется предположение, важное для применения метода конечных элементов к задачам, рассматривающим анизотропные материалы. Главные оси инерции должны быть параллельны координатным осям. Матрицы элемента должны быть составлены относительно главных осей инерции„ которые могут быть различным образом ориентированы относительно глобальной системы координат. Это требование усложняет ввод исходной информации о координатах узлов элемента.
Если координаты узловых точек измеряются относительно глобальной системы координат, то для вычисления координат узлов в локальной системе координат необходимо выполнить преобразование координат. Включение преобразования координат в программу обычно несложная задача, но она требует использования некоторых приемов для того, чтобы указывать. я 82 Глава 8 вычислительной машине момент использования этого преобразог вания.
Эти приемы связаны с вводом дополнительного параметра .для каждого элемента. Запишем формулы преобразования для двумерного случая: (8.58) Здесь х', у' соответствуют локальной системе координат, а х, у— координаты глобальной системы. Эти формулы предполагают совпадение начала координат обеих систем. 8.6. Точечные источники До сих пор мы не рассматривали важные для многих физических задач понятия точечного и линейного источников.
Говорят, что точечный (нли линейный) источник тепла существует, когда генерирование тепла Я происходит внутри очень малого объема илн очень малой ~площади. Физическими и примера~ми линейных источников являются проложенные в земле трубы для подачи горячей воды ~н (илн) водяного дара и проводящий электрический ток провод (х,,у,) в электропроводящей среде. В каждом из этих случаев площадь поперечного сечения трубы или провода, мала по сравнению с,размера~ми окружающей .среды.
В задачах о течении грунтовых ~вод насосы, Фиг. 8лн Точечный источник внутри треугольного еле- выкачивающие;воду из водоносного слоя, менте. также рассматриваются как точечные ис- точники. Точечные и линейные источники достаточно часто встречаются :в окружающей действительности, что оправдывает наше внимание к ним. Мы ограничимся обсуждением источника внутри двумерного элемента, но процедура, которая будет рассмотрена, очень быстро распространяется и на трехмерный элемент. Рассмотрим треугольный элемент на фиг. 8.4 с линейным источником Яе (кВт/м] (тепло притекает и считается положительным), расположенным в точке Хе, Уо. Так как тепловой источник находится в точке, Я теперь не постоянно внутри объема элемента, .а является функцией координат х и у.
Используя единичные имлульсные функции 6(х — хо) и 6(у — уо) 11], можно записать Я =9*6 (х — х,) 6 (р — у,). (8.59) 153 Перенос тепла за счет теплопроеодности и конеекции Интеграл ~(у)т (зл)т может быть теперь записан как т У1 (У)т Яс((т=Яо ~ Ут 6(х — хо) 6Ь вЂ” уо) с(хс(у. (8.60р А Уа Толщина элемента предполагается равной единице. Используя из- вестные свойства импульсных функций, имеем У, (У)т ()с(1 Оо У У к=хо не то (8.61) Соотношение (8.61) устанавливает, что, если точечный (линейный) источник находится внутри элемента, Я распределяется по узлам пропорционально соответствующим величинам Уь У; и Уа, которые вычисляются с помощью координат точечного источника. Так как ХУа=1 в любой точке внутри элемента, мы ие получим величины, большей чем Яо.
Пример 61. Интенсивность источника Яо=52 Вт/см. Источник находится в точке с координатами (5.2) внутри элемента, показанного ниже (этот элемент использовался в задаче 60). Требуется определить распределение яо по узлам элемента. 1ч задаче 61. Глава 8 Значения Ь и с были вычислены Ь,.= — 4, Ь,= Ь = 3, в задаче 60: с,.= — 1, с~ — — — 3, Ь»= 4. Вычислим константы а: а, =Х;1'„— Х,[; =(7) (4) — (6) (О) =28, а~ — — Х»К, — Х;1'» =(6) (3) — (3) (4) = 6, ໠— Х,.У1 — Х;[», = (3) (О) — (7) (3) = — 21.