Метод конечных элементов в задачах теплопроводности (стр. 3 )

• задание размеров и, тем самым, количества элементов;

• нумерацию элементов и узлов, и индексацию последних.

Достаточно подробное рассмотрение этих операций составляет содержание данной главы.

Отметим, что первый этап не связан ни с физической природой краевой задачи, ни с версией – вариационной или взвешенных невязок,− используемой при ее решении. Первый этап – это чисто геометрический аспект метода конечных элементов, имеющий своей целью формирование геометрической части таблицы входных данных, необходимой для машинной реализации МКЭ.

3.1 Типы конечных элементов. Базовый каталог элементов

Согласно основной концепции МКЭ каждый элемент рассматривается независимо от остального их ансамбля (см. п. 2.2). Размерность элемента определяется, очевидно, размерностью аппроксимируемой их совокупностью области определения задачи. Для соответствия элемента физической модели исследуемого объекта одномерный элемент может иметь поперечное сечение, площадь которого не обязательно постоянна по длине элемента, а двумерный элемент – единичную толщину, как это показано на рис. 3.1.

Так как каждый элемент рассматривается независимо от остальных, это дает возможность создать каталог (или архив) разнообразных по геометрии и размерности конечных элементов, обладающих различным порядком симметрии. С помощью такого набора элементов можно собрать любую конструкцию, подобно тому, как это

Рис. 3.1 Двумерная область единичной толщины и одно­мерная

область с произвольным сечением

Согласно основной концепции МКЭ каждый элемент рассматривается независимо от остального их ансамбля (см. § 2.2). Это дает возможность создать каталог (или архив) разнообразных по геометрии и размерности конечных элементов, обладающих различным порядком симметрии. С помощью такого набора элементов можно собрать любую конструкцию, делается с помощью элементной базы детского конструктора. Чем обширнее будет каталог, тем большее разнообразие конструкций может быть реализовано с его помощью.

Краевые задачи теории поля независимо от их физической природы традиционно подразделяют по признаку размерности изучаемого объекта на одно-, двух - и трехмерные. В задачах же теплопроводности при наличии сложного (трехкомпонентного) теплообмена на внешних и внутренних поверхностях объекта такое подразделение с физической точки зрения неприемлемо, так как поверхность (внешняя и внутренняя) самого объекта в этих случаях в теплообмене не участвует – последний осуществляется лишь на внешних границах двумерного объекта, получаемых путем искусственно вводимой его единичной толщины, или на поверхности приписываемого одномерному объекту сечения произвольной формы, как это показано на рис. 3.1. В силу сказанного базовый каталог элементов должен содержать только трехмерные элементы.

В целях упрощения математического описания элементов и сохранения наглядности при декомпозиции объекта на элементы, в качестве таковых следует брать ячейки, образованные координатными поверхностями системы с соответствующим объекту порядком симметрии . Узлами элемента будут точки пересечения координатных линий. Для расширения элементной базы можно использовать и элементы, которые лишь частично образованы координатными поверхностями (трансляцией вдоль одного из ортов системы), либо трансляцией точек в направлениях , не совпадающих с ортами .

Базовый каталог конечных элементов с иллюстрацией их геометрической формы, являющийся рабочим инструментом при декомпозиции области, представлен в Приложении 1. В левом столбце каталога указаны:

• номер элемента , присвоенный ему в каталоге, с указанием орта

трансляции; трансляции;

• порядок симметрии ;

• количество r узлов элемента.

В инженерной практике часто встречаются конструкции, собранные из осесимметричных элементов − цилиндров, конусов и т. д., порядок симметрии которых . Соответственно геометрии объекта задачи теплопроводности с − так называемые осесимметричные задачи, − в общем случае следует разделить на два типа задач:

•1) полностью симметричные − граничные условия не зависят от угловой

координаты θ;координаты ;

•2) ограниченно симметричные.

Ограниченно симметричные задачи с являюется особыми инженерными задачами теплопроводности в силу следующих причин:

•а) объект представляет представляет собой какую-то часть (по углу или углам) осесимметричной фигуры

осесимметричной фигуры;

•б) распределение тепловых − объемных и (или) поверхностных − нагрузок носит локальный характер (граничные условия не являются осесимметричными) в силу особенностей самой конструкции и (или) условий

ее функционирования, (чэто присуще, например, всем космическим объектам [32]);

в) • объект − многосвязная область (например, перфорированные цилиндр, конус, шар и т. д.).

Указанные признаки ограниченной симметрии могут

присутствовать в задаче порознь, в любой комбинации или все одновременно.

В базовом каталоге конечных элементов даны ограниченные по углу (углам) трехмерные элементы , образованные координатными поверхностями системы. Они легко получаются поворотом на угол ∆θ соответствующих плоских фигур (на рисунке заштрихованы). Элементы с номерами 2, 6, 7 лишь частично образованы координатными линиями: путем трансляции треугольника

в направлении или поворотом его на угол вокруг оси . И только третий элемент − тетраэдр – никак не связан с координатной сеткой: он получен трансляцией точек (узлов).

Осесимметричные элементы легко получить поворотом на угол заштрихованных элементообразующих поверхностей (граней).

Включенные в базовый каталог элементы позволяют собрать конструкцию сложной геометрической формы, чем и обусловлено одно из важных преимуществ МКЭ перед другими численными методами.

3.2 Дискретизация области на элементы

В общем случае конструкции различной сложности редко имеюет единый порядок симметрии, – как правило, разные ее части обладает разными величинами . Поэтому в конструкции сначала выделяют части с одинаковым порядком симметрии (необязательно одинаковые, очевидно, по геометрической форме), которые затем разбиваются на конечные элементы, содержащиеся в базовом каталоге с соответствующим порядком симметрии. Процедура представления объемной конструкции совокупностью элементов показпроиллюстрирована на рис. 3.2.

Узел элемента не может располагаться на линии, соединяющей узлы граничащего с ним элемента. Если границами является разные по размерам элементы, то соблюсти указанное требование можно двумя способами: а) при

Рис. 3.2. Разбиение объемной конструкции на элементы

а) с разными (вид сверху); б) с одинаковыми .

Узел элемента не может располагаться на линии, соединяющей узлы граничащего с ним элемента. Если границами является разные по размерам элементы, то соблюсти указанное требование можно двумя способами:

• при сохранении однотипности элементов увеличением их количества;

• б) использованием элемента (элементов) другого типа. Оба способа приведены на рис. 3.3.

Выгодность второго способа очевидна: в первом случае число элементов, аппроксимирующих большеразмерные части конструкции, удвоилось при одновременном возрастании общего количества узлов; во втором случае количество узлов осталось неизменным, а число элементов увеличилось всего на два независимо от того, на сколько элементов была разбита большеразмерная часть конструкции.

Размеры элемента задаются с одновременным учетом следующих условий: на одной и той же поверхности элемента граничные условия должны

Рис. 3.3. Способы разбиения области на элементы разных размеров.

Цифры – это номер элемента по каталогу

Размеры элемента задаются с учетом следующих условий:

• на одной и той же поверхности элемента граничные условия должны

быть физически одинаковыми, т. е. должна соблюдаться равномерность быть физически одинаковыми, – должна соблюдаться равномерность

распределения по поверхности всех видов тепловых нагрузок;

распределения по поверхности всех видов тепловых нагрузок; • объемная нагрузка не может занимать часть элемента;

• материал элемента должен быть одним и тем же по физическим свойствам.свойствам.

Очевидно, что чем меньше размеры элемента, тем с большей точностью выполняется первое требование. Уменьшать размеры элементов желательно в тех частях конструкции, где ожидаются (интуитивно) наиболее резкие изменения искомой функции. Это особенно важно при использовании элементов низкого порядка, обеспечивающих лишь постоянство градиента функции.

Для полностью осесимметричной задачи вследствие отсутствия угловой зависимости по дискретизируется сечение объекта плоскостью, проходящей через ось симметрии. Отсутствие одной из координат понижает на единицу размерность базисных функций, что существенно упрощает процедуры дискретизации области и получения решения (т. к. ). Количество узлов при декомпозиции плоскости, очевидно, много меньше, чем при декомпозиции объема на ограниченные элементы. Поэтому в двумерном случае можно использовать мелкую сетку, т. е. элементы малых размеров.

Использование генетической связи между ограниченно симметричными и осесимметричными элементами позволяет не вводить в базовом каталоге отдельно номера для последних, а ограничиться индексацией номера по "трансляционному" орту , как это показано в каталоге.

Искусство разбиения объекта на элементы зависит от имеющихся навыков. При их отсутствии или недостаточности надеяться на хорошие результаты не приходится. В программах типа “NISA” или “ANSYIS” разбиение на элементы может быть осуществлено в автоматическом режиме, но при этом осуществляется слишком мелкое разбиение, вследствие чего требуется большой объем памяти и длительность времени расчета существенно возрастает. Поэтому при проведении расчетов предварительного характера достаточно проводить дискретизацию конструкции на элементы в ручном режиме.

3.3 Нумерация элементов и узлов

Элементы, на которые разбита конструкция, необходимо индивидуализировать, что проще всего достигается присвоением ему номера . Он никак не связан с номером элемента по каталогу, фиксирующим его геометрию. Нумерация элементов не влияет на вычислительные аспекты МКЭ и поэтому представляет собой простую процедуру, опирающуюся на естественное пожелание удобства при пользовании. Очевидно, что элементы, относящиеся к частям конструкции при укрупненном ее расчленении, должны

иметь последовательную нумерацию. Номер элемента будем заключать в круглые скобки – – во избежание путаницы с номером (по каталогу) и с номерами узлов.

Нумерация узлов существенно влияет на эффективность вычислений. Применение МКЭ к решению дифференциального уравнения приводит к системе алгебраических уравнений (необязательно линейных), большое число коэффициентов в которой равно нулю. Все ненулевые коэффициенты (и некоторые нулевые) в глобальной матрице коэффициентов находятся между двумя линиями, параллельными главной диагонали. Расстояние между главной диагональю и этими линиями называется шириной полосы матрицы. Все коэффициенты вне этой полосы равны нулю, и они не заносятся в память ЭВМ (это одно из преимуществ МКЭ). Уменьшение ширины полосы приводит к уменьшению требуемого объема памяти и к сокращению времени вычислений.

коэффициенты вне этой полосы равны нулю, и они не заносятся в память ЭВМ (одно из достоинств МКЭ). Уменьшение ширины полосы приводит к уменьшению требуемого объема памяти и к сокращению времени вычислений.

В конкретных расчетах структура

матрицы может быть представлена набором целочисленных пар , каждая из которых означает пару переменных (т. е. номера строки и столбца). Полуширина М матрицы определя-а

матрицы может быть представлена набором целочисленных пар , каждая из которых означает пару переменных (т. е. номера строки и столбца). Полуширина М матрицы ется при этом максимумом величины , взятой по всем -элементам матрицы [3]:

. (3.3.1)

При работе с векторными величинами (например, скорость или перемещение в узле), величину нужно умножить на число неизвестных в узле (число компонент векторной величины). В общем случае:

. (3.3.2)

Для скалярной величины, такой как температура, очевидно, .

Объем памяти, необходимой для профильной записи матрицы, определяется формулой:

. (З.3.3)

Правильной нумерацией узлов, очевидно, будет та, которая минимизирует либо полуширину , либо профиль , в зависимости от предполагаемой формы

записи. В большинстве случаев минимизация минимизирует и .

На рис. 3.4 представлены различные варианты нумерации узлов
(и элементов). Сопоставление получаемых показывает предпочтительность последнего варианта, в котором обеспечивается наименьшая из максимальных разница между номерами узлов, принадлежащих одному элементу. Это достигается последовательной нумерацией узлов при движении в направлении наименьшего размера тела. Если от четырехугольных элементов перейти к треугольным, то при правильной нумерации меньшие получаются при проведении левой диагонали, как показано на рис. 3.3. Переход к треугольным элементам лишь удваивает число элементов, не

Рис. 3.4 К вопросу нумерации узлов

На рис. 3.4 представлены различные варианты нумерации узлов
(и элементов). Сопоставление получаемых показывает предпочтительность последнего варианта, в котором обеспечивается наименьшая из максимальных разница между номерами узлов, принадлежащих одному элементу. Это достигается последовательной нумерацией узлов при движении в направлении наименьшего размера тела. Если от четырехугольных элементов перейти к треугольным, то при правильной нумерации меньшие получаются при проведении левой диагонали, как показано на рис. 3.3. Переход к треугольным элеменсказываясь на количестве узлов. Как уже указывалось, нумерация элементов носит произвольный характер, так как формирование глобальных матриц из матриц элементов осуществляется по номерам узлов, а не элементов.

3.4 Индексация узлов и формирование таблицы входных данных

При математическом описании элемента необходимо узлы элемента индивидуализировать узлы элемента, т. е. сделать их отличимыми друг от друга.

С этой целью каждому узлу присваивается индекс в зависимости от количества узлов элемента. Выбор -го узла элемента произволен, индексация остальных узлов выполняется последовательно в направлении против часовой стрелки. Проиндексировав узлы, составляют таблицу, которая ставит в соответствие индексы элемента глобальным номерам его узлов. С помощью этой таблицы впоследствии осуществляется включение матрицы элемента в соответствующие глобальные.

Фиксация положения элемента в пространстве осуществляется заданием координат его проиндексованных и занумерованных узлов.

Процесс дискретизации области определения задачи завершается формированием геометрической части таблицы входных данных.

На рис. 3.54 изображена двумерная (для наглядности) область,

Рис. 3.54 Пример дискре­тизации сечения

осесимметричной детали

разбитая на два элемента. Узлы, выбранные в качестве i-х, обозначены звездочками. В скобках указаны номера элементов и координаты узлов . Геометрическая часть таблицы входных данных с указанием номера элемента по каталогу имеет следующий вид.

Таблица входных данных должна содержать, естественно, сведения о свойствах материала элемента – его механических, теплофизических, электрических и т. п. свойствах, соответствующих физической природе задачи. Эти сведения целесообразно объединить в физический каталог, в котором каждому материалу присвоен номер . Например, девятый номер присвоен ниобию , а третий - нержавеющей стали. Введение в программу постоянно пополняемого физического каталога унифицирует ее и делает компактной таблицу входных данных, поскольку громоздкие, занимающие несколько разрядов данные о свойствах материала элемента (притом для каждого, даже если они одинаковые) заменяются не более чем тремя разрядами номера материала по физическому каталогу.

Задание 3

3.1 Дискретизируйте объект – треугольную прямую призму – на два неодинаковые по размерам элемента. Задайте координаты узлов; составьте таблицу входных данных.

3.2 Дискретизируйте делением по углу и по длине на 4 элемента объект,

Таблица 2

Геометрическая часть таблицы входных данных

Номер элемента	Номер элемента по физическому каталогу	Номер элемента по базовому каталогу	Индекс и номер узла	Координаты узлов
(e)	eфк	ekei	i j k l	Ri Rj Rk Rl	Zi Zj Zk Zl
(1)	9	42	1 3 4 2	0,0 2,0 2,0 0,0	1,0 1,0 3,0 3,0
(2)	3	52	3 5 6 4	2,0 3,0 3,0 2,0	1,0 1,0 3,0 3,0

Геометрическая часть таблицы входных данных

Таблица входных данных должна содержать, естественно, сведения о свойствах материала элемента – его механических, теплофизических, электрических и т. п. свойствах, соответствующих физической природе задачи. Эти сведения целесообразно объединить в физический каталог, в котором каждому материалу присвоен номер . Например, девятый номер присвоен ниобию , а третий - нержавеющей стали. Введение в программу постоянно пополняемого физического каталога унифицирует ее и делает компактной

таблицу входных данных, поскольку громоздкие, занимающие несколько разрядов данные о свойствах материала элемента (притом для каждого, даже если они одинаковые) заменяются не более чем тремя разрядами номера материала по физическому каталогу.

Задание 3

3.1 Дискретизируйте объект – треугольную прямую призму – на
два неодинаковые по размерам элемента. Задайте координаты
узлов; составьте таблицу входных данных.

3.2 Дискретизируйте делением по углу и по длине на 4 элемента
объект,

представляющий собой четверть (по углу ) полого
цилиндра. Задайте координаты узлов; составьте таблицу входных данных.

3.3 Какие свойства материала элемента следует внести в физический каталог в задачах теплопроводности?

ГЛАВА 4. МАТЕМАТИЧЕСКОЕ ОПИСАНИЕ ЭЛЕМЕНТА

Изложенная в п. 2.2 основная концепция метода конечных элементов математически описывалась соотношениями (2.2.3) и (2.2.1):

.

Это означает, что независимо от используемой в дальнейшем версии МКЭ – вариационной или взвешенных невязок, ищется приближенное решение для каждого -го элемента вне связи с решениями для остальных элементов, а затем суммированием этих "элементных" решений находится приближенное решение задачи.

Процедура отыскания на первом этапе заключается в нахождении конкретного вида базисных функций элемента, т. е. в математическом его описании. Подробное рассмотрение этого этапа приводится в данной главе.

4.1 Метод Крамера

Математическое описание элемента можно получить тремя способами. В

4.1 Метод Крамера

Математическое описание элемента можно получить тремя способами. В этом параграфе приведен первый, как наиболее общий.

Для наглядности рассмотрение проведем на конкретном примере двумерного треугольного элемента с тремя узлами (см. рис. 4.13).

Представим приближенное решение для элемента полиномом 1-го ранга с

неизвестными коэффициентами :

, . (4.1.1)

Подставляя в (4.1.1) координаты узлов и получаемые в результате значения функции в каждом из узлов – , соответственно, получим систему уравнений (4.1.2) для определения :

Рис. 4.1

Подставляя в (4.1.1) координаты узлов и получаемые в результате значения функции в каждом из узлов – , соответственно, получим систему уравнений для определения :

=, (4.1.2)

решение которой можно получить методом Крамера.

Определитель этой системы уравнений равен, как известно, 2А – удвоенной площади элемента [11]:

. (4.1.3)

Запишем в развернутом виде решение системы (4.1.2) на языке алгебраических дополнений:

; ;

.

Подставим найденные значения в (4.1.1) и сгруппируем члены, умножаемые на узловые значения функции Uq (q=iI,j,k):

(4.1.4)

Выражения в квадратных скобках зависят от координат узлов элемента и текущих переменных и . Их принято называть функциями формы, базисными или интерполяционными функциями элемента. Представим их в общем виде:

, . (4.1.45)

Элементами матрицы qp служат алгебраические дополнения определи:

(4.1.56)

и являются определителями 2-го порядка.

Более удобной является несколько иная форма записи базисных функций:

; . (4.1.67)

Сопоставляя последнее выражение с (4.1.45) и (4.1.56), видим, что коэффициенты – это столбцы матрицы (4.1.56) при фиксированном соответственно;

при фиксированном это элементы соответствующей строки этой же матрицы. Например, при и получим:

, , ;

при и для базисной функции будем иметь:

.

Раскрывая алгебраические дополнения, найдем конкретные выражения коэффициентов через координаты узлов элемента:

ai = XjJYk –- XkYjJ ; ajJ = XkYi – XiYk ; ak = XiYjJ – XjJYi ;

bi = YjJ – Yk ; bjJ = Yk – Yi ; bk = Yi – YjJ ; (4.1.7)

ci = Xk – XjJ ; cJ = Xi – Xk ; ck = XjJ – Xi .

Переход к другим системам координат осуществляется заменой текущих переменных , : в цилиндрической – на , ; в сферической – на и .

С введением понятия базисной функции аппроксимирующую функцию (4.1.1) (или (4.1.4)) можно представить как явную функцию ее узловых значений :

(4.1.8а))

или в матричной форме: