Дульнев Г.Н., Парфенов В.Г., Сигалов А.В. Применение ЭВМ для решения задач теплообмена (1185899), страница 27
Текст из файла (страница 27)
Отметим, что эти функции существенным образом зависят от вида используемых элементов и способа расположения узлов. При решении двумерных задач Рис. 4.3 применяют элементы различной формы (рис. 4.3, а — г): треугольные, четырехугольные, элементы с криволинейными границами, с узловыми точками, расположенными по границам элементов. Прн решении трехмерных задач элементы могут быть выбраны в виде тетраэдров, косоугольных параллелепипедов, призм с треугольным оса/ Рис. 4.4 нованием (рис. 4.4, а — 6). Ниже ограничимся рассмотрением наиболее простых и часто применяемых линейных треугольных элементов с узлами, расположенными в их вершинах (см.
рис. 4.3, а). Название «линейные» у этих элементов обусловлено тем, что при наличии на границе трех узлов допустимо использовать только функции формы г' (х, у), линейные по координатам х и у: )с(х, у) =а+ 6х+су. (4. 10) Три коэффициента а, Ь н с определяются из трех условий: в двух вершинах значения )с (х, у) равны нулю, а в одной — единице.
Применение более сложных элементов позволяет использовать для 132 аппроксимации температурного поля в элементе полиномы более высоких степеней. Например, для треугольного элемента с шестью узлами (см. Рнс. 4.3, б) нужно применять более сложные посравнению с (4.10) функции формы вида г (х, у)=а+Ьх+су+йху+ехв+(ув. (4.11) Отметим, что при построении функций формы кроме сформулированных выше условий на их значения в узловых точках (единица в одном узле и ноль в остальных) необходимо обеспечивать непрерывность аппроксимации на границах смежных элементов.
Методика получения функций формы для различных видов элементов рассмотрена в 17, 27). Итак, разобьем двумерную область.0 на Лl треугольных элементов и введем М узлов во всех вершинах треугольников (см. рис. 4.2). Причем вершины одних треугольников не должны находиться на сторонах других треугольников, т.
е. каждый элемент должен содержать три узла. Присвоим сквозную нумерацию всем элементам (и = 1, ..., Л() и всем узлам (т = 1, ..., М). Номер узла в общем случае не связан с номером элемента. Вопрос об оптимальном способе нумерации узлов обсудим далее. На рис. 4.5 отмечен некоторый и-й лемент, содержащий верши'- ны с номерами 1, 1, й, которые могут принимать значения от 1 до М. Эти узлы имеют координаты х;, х;, хв по оси х и у;, угь у» по оси у.
Выполненное разбиение удобно описать для дальнейшего использования в программах следующим образом: задать массивы координат всех узлов (х ), и (ум)м=~ в поп м рядке их нумерации; указать для каждого и-го элемента по три номера узлов, лежащих в его вершинах; указать для каждого и-го элемента признак, определяющий при; надлежность его сторон внешней границе области В.
Номера узлов в вершинах и признаки принадлежности сторон границе будем задавать в виде следующей таблицы: Номер элемента Номера узлов Признак принадлежности сторон границе О или 1 или 2 Такую таблицу можно задать в виде матрицы размером Лl Х 4, которую будем называть индексной матрицей. Номера узлов 1, 1, й в индексной матрице будем записывать в порядке, соответствующем обходу элемента против часовой стрелки.
Признак принадлежности сторон границе принимает значение О, если стороны элемента не прилегают к внешней границе; значение 1, если граничной является 133 уу У> к„ хк Рис. 4.5 х к Рис. 4.6 иы>(х, У)=и>Р>ы>(х, У)+и; Р(с>(х„У)+и„Р'">(х у) (4.12) где Р,'">, Р,'"', Р4~" > — линейные функции координат х, у вида 14.9), равные единице в узлах >, 1 или й соответственно и равные нулю в двух других узлах (рис.
4.7). Индекс (л) указывает, что функции формы относятся к л-му элементу. Рис. 4.7 сторона между узлами >' и(; значение 2, если граничными являются две стороны — между узлами > и > и между узлами( и й. Ясно, что других ситуаций можно избежать путем соответствующего выбора узла >', с которого начинается запись номеров в строке индексной матрицы. Предложенный способ полностью описывает выполненное разбиение на элементы и нумерацию узлов и не содержит избыточной информации.
Искомыми величинами в МКЭ являются приближенные значения температуры и в узлах л> = 1, ..., М. Согласно основной идее метода распределение температуры в каждом элементе представляется в виде суммы, в которую входят три функции формы элемента, умноженные на приближенные значения температуры в его трех узловых точках.
Поэтому распределение температуры в л-м треугольном элементе и<'> (х, у) имеет вид (рис. 4.6): Таким образом, для функции формы Р,'"' (х, у) должны выполняться равенства Р,"'(хь У;)=1 Роп(хп у;) =Р<"1(х», у»)=-О. (4.13) Аналогичные соотношения имеют место и для двух других функций Р/"' (х, у), Р»"' (х, у). Используя условия (4.13), можно выразить коэффициенты а'"1, Ь'"1, с1,"1 (т = /, /, /г) через координаты узлов /, /, й и определить функции формы. Приведем выражения для этих коэффициентов. Здесь и далее индекс элемента (п) будем опускать, помня, однако, что функции формы различны для каждого элемента.
Из условий (4.13) следует: а; =(х1 у» — х, у1)/25 а1 =(х„у; — у„х;)/25 ໠—— — (х; у; — х; у»)/25 Ь; = (у1 — у»)/25 Ь; = (у» — у~)/25 Ь» = (у» — у1)/25 с; = (х» — х;)/25 с1 — — (х; — х,)/25 с» =(х; — х»)/25, (4.! 4) где 5 — площадь треугольника, значение которой может быть вычислено через координаты узлов: 5=(х1У» — х»у1+У1х; — у»х~+х„у,— хэу~)/2. (4.!6) Градиент температуры в каждом элементе имеет постоянное значение, и производные по координатам определяются соотношениями: (4.
18) 1ЗЗ ди ди — = Ь» и~ + Ьэ из+ Ь» и,, — = с; и;+ сз и1+ с» и». (4. 16) дх ду При выводе разностных уравнений МКЭ нам понадобится вычис- лять интегралы по площади треугольника и по отдельным сторонам от выражений, содержащих функции формы. Приведем ряд ис- пользуемых в дальнейшем формул. Интеграл по площади элемента от любой его функции формы равен ЦРо'1(х, у)бхбу=5оп/3, т=/, /, /г. (417) зп Как видно из рис.
4.7, этот интеграл равен объему пирамиды с тре- Угольным основанием и высотой, равной 1. Интегралы по стороне Ь;1 вычисляются по формулам: ) Р; (х, у)»)/ = 41/2, ) Р) (х, у)»11 = Т.ы/3, сы сц ) Р;(х, у)Р1(х, у)Й =1.м/6, с,.э где Ец = 1(х» — х;)'+ (у; — уз)»)м» вЂ” длина стороны между уз- лами 1 и /. Данные выражения справедливы для всех функций фоРмы Рм Р;, Р„. $ Л.з. СИСТЕМА УРАВНЕНИЙ МЕТОДА КОНЕЧНЫХ ЭЛЕМЕНТОВ. ЛОКАЛЬНАЯ И ГЛОБАЛЬНАЯ МАТРИЦЫ Перейдем к следующему этапу — составлению системы разностных уравнений !)4КЭ.
Система уравнений для определения температур в узлах (и„) ! составляется на основании условий минимум ма функционала (4.8). Этот функционал можно представить в виде суммы интегралов по всем элементам, на которые разбита область е): — Т(и) (4. 19) л= ! У(л) = ~ [Л ~ — ) +- Л( — ) — 2() и~бх((у+ з (л) + ~ (аи' — 2()ли)(1!. (4. 20) с(л) Второй интеграл в (4.20) вычисляется лишь для граничных элементов вдоль тех их сторон, которые прилегают к внешней границе !. области В. Условия минимума функционала можно с учетом (4.19) переписать в виде д (~~ у(и) ~ у д!л 0 т 1 ((4 (42!) дил, ( и ( / лле дл!и л=! л=! (36 Структура системы разностных уравнений.
Рассмотрим подробнее структуру системы разностных уравнений (4.21). Возьмем т-е уравнение, получающееся при дифференцировании функционала (4.19) по значению температуры в т-м узле. Напомним, что распределение температуры и(л) (х, у) в любом элементе зависит только от температур и;, из, ид в узлах этого элемента. Соответственно и значение функционала!(л) зависит только от этих температур.
Поэтому в сумме (4.!9) от и будут зависеть только !(л) тех элементов, которые включают т-й узел. Это обстоятельство позволяет формировать систему (4.21) одним из двух способов. Первый способ основан на переборе узлов. Он заключается в выделении какого-то т-го узла, определении элементов, содержащих этот узел, и суммировании частных производных от функционалов этих элементов по температуре выделенного узла. Таким образом, при получении системы (4.21) последовательно формируется первое, второе, ..., ()4-е уравнения. Второй способ основан на переборе элементов.
В этом случае выделяется какой-то л-й элемент и три производные от функциона- д l ди >» — ~ — ) = 2 (Ь( и(+ Ь; и> + Ь„и«) Ь(, д«п дх ) д I ди>« — ~ — ) =2(с; и(+сти;+с„и») с>, " ~ь1 — (2(>„и) = — 2(>, Р((х, у), д д>п — (2(),и) =2(ЬР,(х, у), д д«п — (аи') =- 2а (Р( и; + Р; и;+ Р» и») Р(.
д И( (4. 22) Будем считать, что разбиение выполнено таким образом, чтобы величины Л, д„, ()„а можно было бы считать постоянными в пределах каждого элемента. Если эти величины являются кусочно-непрерывными функциями координат, то разбиение проводят так, чтобы границы элементов совпадали с линиями разрыва. Если величины являются непрерывными, но резко изменяющимися функциями координат, то нужно построить в области их сильного изменения более густую сетку.
Вычислим сначала интеграл по площади элемента, используя выражение (4.17) для интеграла от функций формы: — ~Л ( — ) + Л( — ) — 2(),и~([х([у= з(«) ) 2 [Л (Ь( и; -(- Ь> и + Ь„ид) Ь(+ Л (с> и(+ с; ит+ с» и») с(— з(л) — д, Р((х, у)[([х((у=25("> [Л(Ь»и>+Ь>Ь;и>+Ь(Ь„и»)+ + Л (сз и(+ с( сд и>+ с( с„и„) — (),/3[. (4.23) 137 ла этого элемента по его узлам заносятся в левые части соответствующих трех уравнений системы (4.21), т.
е. производится формирование сразу всей системы, которое реализуется путем постепенного «наращивания» левых частей уравнений. Этот способ получил более широкое распространение и рассматривается ниже. Очевидно, что для его реализации понадобятся выражения для производных от функционала и-го элемента !("> по температурам его узлов и(, иь и„.,([ля их получения поменяем порядок дифференцирования и интегрирования, т. е. сначала проведем дифференцирование подынтегрального выражения в (4.20), а затем вычислим интеграл по элементу от получившейся функции. Прн реализации этой процедуры будем использовать выражение (4.12) для распределения и("> (х, у), причем индекс (и) для сокращения записей опустим. Учитывая выражения (4.16) для производных, имеем Если стороны элемента принадлежат границе области, то в функционале следует учесть интеграл по этим сторонам.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
