Методом конечных элементов в mathcad. Решение двумерной задачи теплопроводности методом конечных элементов в mathcad Функционал решения задачи теплопроводности

Раздел 7. Метод конечных элементов (МКЭ).

7.1. ЗАДАЧА. Выполнить расчет сил, действующих в узлах элемента (рис. 7.1) .

Рис. 7.1. Шарнирно опертая балка C
p - распределенная поперечная нагрузка, E 1 - модуль упругости балки, A 1 - постоянное сечение балки, L - длина балки, x i , y i - координаты узловых точек балки, u i , v i - перемещения узловых точек балки, U i , V i - перемещения балки, i, n - узловые точки балки

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

На рис. 7.2 изображена двумерная конструкция, состоящая из отдельных частей, соединенных между собой в точках, пронумерованных от 1 до n. Соединения в узлах предполагаются шарнирными.

Рис. 7.2. Типичная конструкция, составленная из отдельных элементов
Сначала допустим, что в результате расчета или на основе экспериментальных данных достоверно известны характеристики каждого элемента. Силы, возникающие в узлах 1 - 3 элемента a , однозначно определяются перемещениями этих узлов, действующей на элемент распределенной нагрузкой p и его начальной деформацией. Начальная деформация может быть обусловлена температурным воздействием, усадкой или несовершенством сборки. Силы и соответствующие им перемещения определяются компонентами U, V и u, v в какой-либо системе координат.

Записывая силы, действующие во всех (в трех для рассматриваемого случая узлах элемента a, в виде матрицы), получим

а для соответствующих перемещений узлов

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

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

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

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

Матрица [k] a называется матрицей жесткости элемента, а [S] a - матрицей напряжения элемента.

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

где F i и б i имеют одинаковое число компонент или степеней свободы.

Ясно, что матрицы жесткости элемента всегда будут квадратными вида

где k ii и т.д. - также квадратные подматрицы размерности l x l , а l - число компонент силы в рассматриваемых узлах.

В качестве примера рассмотрим двумерную задачу о шарнирно опертой C балке постоянного сечения A с модулем упругости E (рис. 7.1). Балка нагружена равномерно распределенной поперечной нагрузкой p и подвержена однородной температурной деформации

г 0 =aT

Если концы балки имеют координаты x i , y i и x n , y n то ее длина может быть вычислена как

А ее угол наклона к горизонтальной оси

В каждой узловой точке необходимо рассмотреть только по две компоненты силы и перемещения.

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

Элементы этой матрицы равны соответствующим компонентам реакций опор балки, то есть pL/2 . Для компенсации температурного расширения г 0 нужно приложить осевую силу EaTA , компоненты которой

Наконец, перемещения узловых точек элемента

вызовут его удлинение (u n -u i)cosa+(v n -v i)sina . Величина удлинения, умноженная на EA/L , даст осевую силу, компоненты которой можно найти, подставив величину этой силы вместо в предыдущее выражение. Стандартная форма имеет вид

Итак, для рассматриваемого простейшего случая определены все слагаемые основного уравнения (7.1). Нетрудно записать в форме (7.2) и напряжения в любом поперечном сечении элемента. Если, например, ограничиться рассмотрением среднего сечения балки C, напряжения, возникающие в результате осевого растяжения и изгиба элемента, можно записать в виде

Где d - половина высоты сечения, а I - момент инерции. В это выражение входят все слагаемые формулы (7.2).

Для более сложных элементов требуются более тонкие приемы расчета, но все равно результаты имеют такую же форму.

Типичная процедура решения задачи МКЭ по методу перемещений для плоского треугольного элемента при упругом поведении материала может быть представлена следующими этапами :

1) Расчетная область разделяется воображаемыми линиями на конечные элементы (КЭ) прос ой формы, например, треугольники.

2) Производится нумерация КЭ и узлов, КЭ связываются между собой в узловых точках, определяются неизвестные в узлах и степени свободы узлов; в качестве неизвестных выбираются перемещения узлов.

3) Выбираются функции (линейные полиномы), аппроксимирующие перемещения в каждом КЭ, которые выражаются через узловые перемещения.

Рассматривается КЭ e , имеющий три узла i, j, m . На рис. 7.4 показан типичный треугольный элемент с узлами i, j, m , пронумерованными против часовой стрелки. Перемещения каждого узла имеют две компоненты

А шесть компонент перемещений элемента образуют вектор

Перемещения внутри элемента должны однозначно определяться этими шестью величинами.

Рис. 7.4. Элемент сплошной среды для расчета плоского напряженного или плоского деформированного состояния
Простейшим представлением являются линейные полиномы

Значения шести постоянных a i легко найти из двух систем, состоящих из трех уравнений, которые получаются в результате подстановки в (7.3) узловых координат и приравнивания значения перемещений соответствующим перемещениям узловых точек. Записав, например,

выражают a 1 , a 2 , a 3 через величины узловых перемещений u i , u j , u m и окончательно

остальные коэффициенты получаются циклической перестановкой индексов i, j, m , а величина определяется соотношением

Аналогично можно представить перемещение v в вертикальном направлении

Соотношения (7.5а) и (7.6) в стандартной форме определяют перемещения любой точки внутри элемента

Где I - единичная матрица размерности 2x2, - координатные функции, которые называются функциями формы

4) Через узловые перемещения выражаются также деформации и напряжения.

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

Используя равенства (7.7) или (7.5а) и (7.6), имеют

что явным образом определяет матрицу [B] .

Матрица [B] не зависит от координат точки внутри элемента, и, следовательно, деформации в нем постоянны.

Матрица упругости [D] , входящая в соотношение, которое в рассматриваемом случае имеет вид

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

Плоское напряженное состояние в изотропном материале. Для плоского напряженного состояния изотропного материала имеем по определению

Разрешая эти соотношения относительно напряжений, получают матрицу [D] в виде

где E - модуль упругости, a v - коэффициент Пуассона. Плоское деформированное состояние в изотропном материале. В этом случае, кроме трех компонент напряжения, существует нормальное напряжение Q z . Для частного случая изотропного теплового расширения

и, кроме того,

Исключая Q z , определяют три остальные компоненты напряжения. Полагая начальную деформацию равной нулю и сравнивая с соотношением (7.10) получают матрицу [D] в виде

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

Для КЭ e вектор-столбец усилий имеет вид

Где - вектор-столбец перемещений узлов данного КЭ, индексы р и о относятся к распределенным и начальным нагрузкам соответственно, - матрица жесткости КЭ

Матрица жесткости элемента ijm определяется с помощью общего соотношения

Где t - толщина элемента, а интегрирование производится по площади треугольника. Если предположить, что толщина элемента постоянна, что тем ближе к истине, чем меньше размеры элемента, то, поскольку ни одна из матриц не содержит x или y , имеют простое выражение

Где - площадь треугольника [введенная соотношением (7.5в)]. Такая форма записи позволяет вычислить матрицу с помощью ЭВМ. Матрицу [B] , определенную соотношением (7.9), можно записать в виде

Матрица жесткости элемента может быть записана в виде

где подматрицы размерности 2x2 строятся следующим образом:

6) Составляется ансамбль КЭ и формируется глобальная матрица жесткости [K] всей расчетной схемы

7) Составляется

И решается система линейных алгебраических уравнений

Считаем узловые силы, обусловленные начальной деформацией и напряжениями, нулевыми.

В общем случае плоского напряженного или деформированного состояния на каждый элемент единичной площади в плоскости x, y действуют распределенные объемные силы

В направлениях соответствующих осей.

Вклад этих сил в узловые силы определяется выражением

или на основании (7.7)

при условии, что объемные силы X и Y постоянны. Так как N i не является постоянной, должно быть выполнено интегрирование.

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

И, используя (7.8), получают

или

Для всякого элемента

Это означает, что все объемные силы, действующие в направлениях x и y, распределены между тремя узлами поровну.

8) По найденным значениям перемещений в каждом элементе определяются в соответствии с постановкой задачи деформации, а затем напряжения.

Ниже реализована в математическом пакете Mathcad рассмотренная процедура решения для тестовой задачи растяжения пластины из двух КЭ. Узловые нагрузки распределяются по всем узлам пропорционально числу КЭ в узле.

Я уже давно использую Mathcad в расчетах, наверное, лет 15, если не больше. Работал в разных версиях Mathcad (2, 5, 6, 7, 8 … 15), то есть практически во всех версиях. Вначале с трудом осваивал, затем с восхищением использовал в расчетах (учебных и научно-исследовательских).

Семь лет назад стал бета-тестером фирмы Mathsoft (принимал участие в тестировании новых версий Mathcad ). Многократно выступал с предложениями по усовершенствованию пакета Mathcad , превращению его в мощное средство программирования и решения сложных расчетных задач. Получал от Mathsoft благодарности и стандартные ответы типаNous hardly thought on those questions (мы напряженно думаем над этими вопросами).

Фирма Mathsoft хоть что-то делала для совершенствования математического аппарата Mathcad . Несколько лет назад фирму Mathsoft купила фирма PTC , основной продукт которой пакет ProEngeneer – ProMechanica . Для них Mathcad побочный продукт. Специалисты, разрабатывавшие Mathcad , ушли. На мой взгляд, изменилась стратегия. Разработанный РТС новый вариант (не версия, а принципиально новый пакет) Mathcad Prime похож на мощный калькулятор с принципиально другим интерфейсом и, на мой взгляд, неудобный в работе (при решении достаточно сложных задач со сравнительно большой программой). К тому же Mathcad Prime пока ещё несовместим со старыми версиями (даже с mathcad 15) По просьбам пользователей РТС вернулась (на время) к старому направлению и выпустила Mathcad 15, являющийся математически 100% копией Mathcad 14. Всё. Приехали.

В каком же состоянии остался, брошенный хозяевами Mathcad ?

Отличный математический пакет для решения задач средней сложности и практически не пригоден для решения сложных задач, типа метода конечных элементов.

Основное препятствие на пути решения сложных универсальных программ отсутствие возможности выбора вариантов расчета . В Mathcad очень сложно обойти ненужный в данном варианте расчета оператор. Ну хотя бы позволили переход на метку, презираемый программистами, но нужный для расчетчиков. (Обойти отдельный оператор можно с помощью оператора ON ERROR).

Что касается метода конечных элементов, то тут я полностью согласен с фирмой РТС. Использование Mathcad для решения задач методом конечных элементов это, извините, “шизо”. Для этого существует много различных вычислительных комплексов, например, ANSYS или тот же ProEngeneer . Я использую Mathcad для решения задач МКЭ только для изучения алгоритма МКЭ. Фактически, выложенная для скачивания курсовая работа представляет собой действующую модель вычислительного комплекса в разрезе. В курсовой работе наглядно виден алгоритм МКЭ (в виде работающих формул). Изучение и запоминание этого алгоритма и есть цель курсовой работы по МКЭ, а вовсе не получение близких к достоверным результатов расчета.

Математическое ядро Mathcad несовершенно (и не совершенствуется вовсе) .
В результате при расчете курсовой работы, например, при корректировке размеров поперечных сечений, при изменении одного из размеров вдруг перестает решаться система уравнений (деление на ноль) или первая собственная частота вдруг оказывается комплексным числом. Небольшая корректировка того же размера и решение появляется вновь . При решении систем уравнений результат расчета иногда резко меняется при изменении значений начальных приближений, или при выборе другого метода расчета (в контекстном меню). Иногда удается повысить точность расчета, выбрав в менюИнструменты (Tools ) -Опции (Worksheet Options ) TOL = 10 -8 вместо TOL = 10 -3 , но далеко не всегда это помогает правильно решить задачу

Министерство образования и науки Российской Федерации

Государственное образовательное учреждение

высшего профессионального образования

«Тихоокеанский государственный университет»

Методические указания и контрольные задания к выполнению

лабораторной работы по курсу «Аналитические и численные методы решения уравнений математической физики» для студентов, обучающихся в магистратуре Хабаровск Издательство ТОГУ 2011 УДК 539.3/6. (076.5) Решение двумерной задачи теплопроводности методом конечных элементов в MATHCAD: методические указания и контрольные задания к выполнению лабораторной работы по курсу «Аналитические и численные методы решения уравнений математической физики» для студентов, обучающихся в магистратуре / сост. Л. М. Иванников. – Хабаровск: Изд-во Тихоокеан. гос. ун-та, 2011.

Методические указания составлены на кафедре «Механика деформируемого твердого тела». Включают содержание лабораторной работы и рекомендации к изучению разделов курса «Аналитические и численные методы решения уравнений математической физики», необходимых к ее выполнению, список рекомендуемой литературы и задачи для лабораторной работы.

Печатается в соответствии с решениями кафедры «Механика деформируемого твердого тела» и методического совета института строительства и архитектуры.

© Тихоокеанский государственный университет,

ОБЩИЕ ПОЛОЖЕНИЯ

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

УРАВНЕНИЕ ПЕРЕНОСА ТЕПЛА

Уравнение плоской задачи теплопроводности имеет вид 2T (x, y) 2T (x, y) Ky Q(x, y) 0, (1) Kx x 2 y кВт где K x, K y – коэффициенты теплопроводности в направлении осей x, y, ;

(м К) T (x, y) – искомая функция температуры; Q(x, y) – источник тепла внутри тела, кВт. Q(x, y) 0, если тепло подводится к телу.

м Граничные условия ставятся двух типов :

T TГ (Г), 1. (2) если температура T известна на некоторой части границы Г, где TГ (Г) – известная температура в точках границы, зависящая от координат точек поверхности s на границе Г;

T (x, y) T (x, y) l Ky m h(T T) q(x, y) 0, 2. (3) Kx x y если на части поверхности Г 1 происходит конвективный теплообмен, характеризуемый величиной h(T T), или задан поток тепла q(x, y) на части поверхности Г 2, причем Г Г1 Г 2. Обозначения в (2) и (3): h – коэффициент кВт теплообмена, ; T (x, y) – неизвестная температура на границе, К; T – (м2 К) известная температура окружающей среды, К; l, m – направляющие косинукВт сы; q(x, y) – известный поток тепла, считается положительным, если м тепло теряется телом. Поток тепла и конвективная теплоотдача на одном и том же участке не могут действовать одновременно.

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

dT 0, dn где n – внешняя нормаль к границе рассматриваемой области.

ФУНКЦИОНАЛ РЕШЕНИЯ ЗАДАЧИ ТЕПЛОПРОВОДНОСТИ

Решение уравнения (1) по области s с граничными условиями (2) и (3) на Г эквивалентно отысканию минимума функционала При решении задачи МКЭ область s разбивается на n подобластей (конечных элементов), которые обычно принимаются в форме треугольников (рис. 1). Далее все формулы приводятся для треугольных КЭ. Функционал записывается как сумма вкладов всех конечных элементов по области. Тогда (4) примет вид проводности.

Или Представим температуру, изменяющуюся в пределах КЭ, через узловые значения:

где [ N (e) ] – матрица функций формы КЭ, учитывающая распределение температуры в пределах КЭ.

Тогда где [ B (e) ] – матрица градиентов функций формы КЭ.

Для каждого КЭ теперь можно записать вклад каждого КЭ в выражение для функционала (4):

Минимум функционала (4) требует выполнения следующего условия:

Для отдельного КЭ получим где матрица теплопроводности КЭ имеет вид а вектор внешнего воздействия будет Для всей рассматриваемой области получим или где Уравнение (6) является основным уравнением для решения задачи теплопроводности методом конечных элементов.

ДВУМЕРНЫЙ СИМПЛЕКСЭЛЕМЕНТ

Для решения плоской задачи теплопроводности используется треугольный КЭ с прямолинейными сторонами (см. рис. 1). Нумерация узлов проводится против часовой стрелки, начиная с некоторого узла, обозначаемого единицей.

Нумерация сторон КЭ приведена на рис. 1.

Узловые значения температуры обозначаются T1, T2, T3. Температура в точке КЭ с координатами x, y определяется по формуле Ниже приводятся функции формы, применяемые для этого КЭ.

Площадь КЭ вычисляется по известной формуле Коэффициенты, входящие в функции формы, зависят от координат узлов, они приведены ниже:

ПРИМЕНЕНИЕ ЧЕТЫРЕХУГОЛЬНЫХ КЭ ДЛЯ ГЕНЕРАЦИИ СЕТКИ

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

На каждой стороне КЭ вводится по три узла.

На рис. 2 показаны местные относительные координатные оси, в которых узел 7 (1, 1). Нумерация узлов такого КЭ, начиная с узла 1, проводится против часовой стрелки. Узлы 2, 4, 6, 8 могут располагаться в произвольной точке соответствующей стороны, что позволяет в дальнейшем строить более густую сетку вблизи точечных воздействий. В дальнейшем каждая сторона такого КЭ разбивается на заданное число участков. Нумерация узлов проводится следующим образом: по вертикали от узла с координатами (1, 1) вниз по оси и слева направо по оси. Таким образом, крупные элементы делятся на более мелкие, которые в свою очередь меньшей по длине диагональю разбиваются на треугольные КЭ. Треугольные участки зоны также представляются в виде четырехугольных квадратичных элементов (рис. 3).

Рис. 3. Представление треугольной области в виде четырехугольного квадратичного элемента

МАТРИЦА ТЕПЛОПРОВОДНОСТИ КЭ

Для треугольного КЭ матрица теплопроводности имеет вид где L1 2, L2 3, L31 – длины соответствующих сторон КЭ. Последние три члена учитывают конвективный теплообмен по каждой стороне КЭ. Так как КЭ входит составной частью в рассматриваемую область, то конвективный теплообмен обычно происходит по одной или двум сторонам КЭ.

ВЕКТОР ВНЕШНИХ ВОЗДЕЙСТВИЙ НА КЭ

Внешними (известными) воздействиями являются:

1. Источник тепла внутри КЭ постоянной интенсивности Q (e).

2. Приток тепла за счет теплового потока q (e).

3. Конвективный теплообмен не более чем по двум сторонам КЭ с коэффициентом теплообмена h (e).

4. Точечный источник тепла Q * (X 0, Y0), находящийся внутри КЭ.

Вектор внешних воздействий на КЭ имеет вид

ГРАДИЕНТЫ ТЕМПЕРАТУР И СРЕДНЯЯ ТЕМПЕРАТУРА ПО КЭ

Градиенты температур и средняя температура по КЭ вычисляются по следующим формулам:

ПОРЯДОК РЕШЕНИЯ ЗАДАЧИ ТЕПЛОПРОВОДНОСТИ

НАНЕСЕНИЕ СЕТКИ УЗЛОВ НА РАССМАТРИВАЕМУЮ ОБЛАСТЬ

Область решения задачи помещается в систему глобальных координат X, Y. Рассматриваемая область должна быть покрыта сеткой узлов. Чем меньше будет ячейка сетки, тем точнее будет решение задачи. Нанесение сетки проводится согласно в 2 этапа.

I этап. Рассматриваемая область разбивается на ряд прямоугольных и треугольных зон (четырехугольные квадратичные элементы). Зоны нумеруются в произвольном порядке. Для каждой такой зоны задаются 8 узловых точек (по три на каждой стороне, включая угловые точки). Для треугольной зоны одна из сторон соответствует двум сторонам прямоугольника (5 точек).

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

Составляются следующие таблицы исходных данных:

а) Табл. 1 соединения зон, определяющая, какие стороны зон контактируют между собой.

Соединение зон в рассматриваемой области. Таблица 1.

В приведенной табл. 1 показано, что зона 1 контактирует только с зоной по первой стороне, зона 2 контактирует с зоной 1 по первой стороне и с зоной 3 по четвертой стороне. Зона 3 контактирует только с зоной 2 по второй стороне (рис. 4). Нумерация сторон зависит от ориентации местных осей в относительных координатах, которые показаны на рисунке жирными цифрами. На рис. 4 показано направление нумерации узлов зон от начального узла Н.

Рис. 4. Формирование таблицы соединения зон б). Табл. 2 координат узлов, нанесенных на границы зон, в принятой глобальной системе координат.

Координаты узлов на границах зон Таблица 2.

в). Табл. 3, в которой указывается число полос по вертикали и горизонтали, на которые разбивается каждая зона для получения сетки с ячейками меньших размеров.

Формирование сетки с меньшими по размеру ячейками Таблица 3.

Зона 1 разбивается на пять полос по высоте и шесть полос по ширине.

г). Табл. 4, в которой для каждой зоны указываются ранее нанесенные узлы.

Номера узлов предварительной сетки для каждой зоны Таблица 4.

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

II этап. Далее в Mathcad реализована программа “grid” , в которой задается число полос по высоте и ширине для каждой зоны, позволяющее разбить каждую зону на прямоугольники гораздо меньших размеров. Затем каждый из этих малых прямоугольников меньшей по длине диагональю делится на два треугольника и вся рассматриваемая область покрывается сеткой с треугольной ячейкой.

В результате работы этой программы выдаются следующие данные:

a). Число треугольных КЭ (Kol_Elem).

б) Следующие табл. 5, 6, 7.

Нумерация узлов сетки по сторонам зон Таблица Таблица выдается в форме матрицы размером (число полос зоны по высоте число полос зоны по ширине) для каждой зоны, что упрощает построение сетки.

Приведенная матрица показывает, что в зоне 3 на стороне 1 располагаются узлы 23, 24, 25, 26; на стороне 2 располагаются узлы 26, 22, 1; на стороне 3 – узлы 1, 16, 13, 10; на стороне 4 узлы 10, 19, 23. Обход зоны против часовой стрелки. Эта нумерация показана в приведенном ниже примере.

Расположение КЭ и принадлежность узлов КЭ треугольной сетке Таблица Также могут быть выведены таблицы, связывающие номер зоны, номер КЭ и координаты узлов КЭ.

На схему рассматриваемой области вручную наносится сетка с нумерацией КЭ и их узлов.

ФОРМИРОВАНИЕ ВЕКТОРА ВНЕШНИХ ВОЗДЕЙСТВИЙ

На основании построенной сетки для рассматриваемой области отмечаются:

а) Номера сторон, по которым происходит конвективный обмен тепла.

б) Номера узлов, в которых температура задана.

в) Номера КЭ, в которых на их сторонах, узлах или внутри располагаются сосредоточенные тепловые источники.

Составляются следующие табл. 8, 9, 10.

Стороны области с конвективной теплоотдачей Таблица Предполагается, что конвективный теплообмен возможен только по двум сторонам КЭ из трех.

Таблица точечных источников тепла Таблица Таблица величин температуры в узлах КЭ.

Таблица градиентов температур Gradx, Grady по осям Х и Y соответственно.

Таблица средней температуры Тsred по каждому КЭ.

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

ПРИМЕР ВЫПОЛНЕНИЯ ЛАБОРАТОРНОЙ РАБОТЫ

В теплопроводящей среде проходят 4 кабеля, как показано на рис. 5. Среда имеет коэффициенты теплопроводности K x K y 10. Коэффициент тепсм К лообмена на поверхности среды h 5. По боковым сторонам рассматрисм 2 К ваемая среда ограничена толстым слоем изоляции. Температура воздуха на поверхности среды T 30 0 C. Температура нижнего слоя среды T 20 0 C.

Мощность излучения тепла каждым кабелем составляет Q 200 Вт.

Требуется:

УКАЗАНИЯ:

а) при выполнении лабораторной работы учесть симметрию области и симметрию температурного воздействия;

б) разбить рассчитываемую часть области на три или четыре зоны;

в) каждую зону разбивать от трех до пяти полос по высоте и ширине для упрощения нанесения сетки на область.

РЕШЕНИЕ ЗАДАЧИ

Учитывая симметрию рассматриваемой области, в расчете будем учитывать только половину этой области (рис. 6).

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

Рис. 7. Предварительная разбивка области на зоны Для более точного решения задачи необходимо узлы на границе зон располагать ближе к точечным источникам тепла.

Составляются исходные данные для назначенных зон и узлов (табл. 1, 2, 3, 4). Программа расчета выдает табл. 5, 6, 7, представляющие полную информацию о треугольной сетке, нанесенной на область, используемую в дальнейшем расчете. По этим таблицам на листе строится сетка (рис. 8).

Рис. 8. Треугольная сетка, нанесенная на область По полученной сетке проводится учет внешнего температурного воздействия и составляются табл. 8, 9, 10. После чего в табличной форме выводятся результаты решения задачи и их графическое представление на рис. 9 и

1. Р Е З У Л Ь Т А Т Ы РЕШЕНИЯ ПО СОЗДАНИЮ СЕТКИ КЭ

УЗЛЫ СЕТКИ ПО ГРАНИЦАМ ЗОН

ТАБЛИЦА КЭ

КООРДИНАТЫ УЗЛОВ КЭ

ФОРМИРОВАНИЕ ВЕКТОРА

ВНЕШНИХ ВОЗДЕЙСТВИЙ

stor

2. Р Е З У Л Ь Т А Т Ы РЕШЕНИЯ ЗАДАЧИ

Градиенты температур и средняя температура по КЭ пературы по области

ВАРИАНТЫ ЗАДАНИЙ ДЛЯ ЛАБОРАТОРНОЙ РАБОТЫ

В теплопроводящей среде, как показано на схеме, проходят кабели, излучающие тепло. Среда имеет коэффициенты теплопроводности K x и K y. Коэффициент теплообмена на поверхности среды h. На некоторых участках рассматриваемая среда ограничена толстым слоем изоляции. Температура воздуха на отдельных участках среды, где происходит конвективный теплообмен, Т. На некоторых участках среды задана температура Т.

Мощность излучения тепла каждым кабелем составляет Q.

Требуется, используя исходные данные для своего варианта и схему задания (табл. 11, рис. 11):

1. Определить распределение температуры в заданной области.

2. Определить градиенты температур и среднюю температуру по области.

3. Построить графики изменения полученных величин.

ИСХОДНЫЕ ДАННЫЕ

Исходные данные к лабораторной работе по вариантам Таблица Ноh, анта Рис. 11. Схемы вариантов задний для лабораторной работы

КОНТРОЛЬНЫЕ ВОПРОСЫ

Запишите уравнение теплопроводности для двумерной задачи.

Запишите граничные условия для двумерной задачи теплопроводности.

Запишите полный функционал решения задачи теплопроводности.

Выведите основное уравнение для решения двумерной задачи теплопроводности методом конечных элементов.

5. Какие конечные элементы используются для решения двумерной задачи теплопроводности?

6. Как определяются функции формы для двумерного симплексэлемента?

7. С какой целью используются четырехугольные квадратичные элементы?

8. Как выбирается система местных координат и проводится нумерация сторон четырехугольного квадратичного элемента?

9. Запишите матрицу теплопроводности для треугольного КЭ.

10. Как формируется матрица теплопроводности для рассматриваемой области?

11. Как формируется вектор внешних тепловых воздействий для КЭ?

12. Как формируется вектор внешних воздействий для рассматриваемой области?

13. Как определяются градиенты температур и средняя температура по КЭ?

14. Как проводится нанесение сетки на рассматриваемую область?

15. Какие исходные данные необходимо подготовить для нанесения сетки?

16. Какие выходные данные используются для построения сетки и как она наносится на область?

17. Какие данные необходимо внести для формирования вектора внешних тепловых воздействий?

18. Как учесть знак величины точечного источника тепла? Притока тепла?

19. Какие выходные данные получаются в результате решения задачи теплопроводности?

1. Зенкевич О. Метод конечных элементов в технике / О. Зенкевич. – М. :

Мир, 1975. – 452 с.

2. Сегерлинд Л. Применение метода конечных элементов/ Л. Сегерлинд. – М.:

Мир, 1979. – 392 с.

ОБЩИЕ ПОЛОЖЕНИЯ ………………………………………………. Уравнение переноса тепла……………….…………………………… … Функционал решения задачи теплопроводности……………………..... Двумерный симплекс элемент…………………..………………..…. Применение четырехугольных КЭ для генерации сетки....……….… Матрица теплопроводности КЭ…………………….………………... Вектор внешних воздействий на КЭ….…………………………… …... Градиенты температур и средняя температура по КЭ…………… …… Порядок решения задачи теплопроводности в Mathcad 14…..….… ПРИМЕР ВЫПОЛНЕНИЯ ЛАБОРАТОРНОЙ РАБОТЫ…….… ….. РЕШЕНИЕ ЗАДАЧИ……………………………………………………. Распечатка решения задачи……………………………………… …..... ВАРИАНТЫ ЗАДАНИЙ ДЛЯ ЛАБОРАТОРНОЙ РАБОТЫ…..... Контрольные вопросы…………………………………………… ……. Библиографический список……………………………………… ……

РЕШЕНИЕ ДВУМЕРНОЙ ЗАДАЧИ ТЕПЛОПРОВОДНОСТИ

МЕТОДОМ КОНЕЧНЫХ ЭЛЕМЕНТОВ В MATHCAD

Методические указания и контрольные задания к выполнению лабораторной работы по курсу «Аналитические и численные методы решения уравнений математической физики» для студентов, обучающихся в магистратуре.

Главный редактор А. А. Суевалова Редактор Т. Ф. Шейкина Оператор компьютерной верстки Л. М. Иванников Подписано в печать Бумага писчая. Гарнитура «Таймс». Печать цифровая.

Усл. печ. л. Тираж 50 экз. Заказ Издательство Тихоокеанского государственного университета.

680035, Хабаровск, ул. Тихоокеанская, 136.

Отдел оперативной полиграфии издательства Тихоокеанского государственного университета. 680035, Хабаровск, ул. Тихоокеанская, 136.

Похожие работы:

«В.Б. Пономарев А.Е. Замураев АСПИРАЦИЯ И ОЧИСТКА ПРОМЫШЛЕННЫХ ВЫБРОСОВ И СБРОСОВ Федеральное агентство по образованию ГОУ ВПО Уральский государственный технический университет–УПИ В.Б. Пономарев А.Е. Замураев АСПИРАЦИЯ И ОЧИСТКА ПРОМЫШЛЕННЫХ ВЫБРОСОВ И СБРОСОВ МЕТОДИЧЕСКИЕ УКАЗАНИЯ ПО КУРСУ МАШИНЫ И АГРЕГАТЫ ПРЕДПРИЯТИЙ СТРОИТЕЛЬНЫХ МАТЕРИАЛОВ Научный редактор – проф., канд. техн. наук В.Я.Дзюзер Екатеринбург УДК 666.9.001.575 (042.4) ББК 35.41в П Рецензенты: Пономарев В.Б. П56 Аспирация и...»

«Министерство образования и науки Российской Федерации Федеральное государственное бюджетное образовательное учреждение высшего профессионального образования Тихоокеанский государственный университет ПРОИЗВОДСТВО КАМЕННЫХ РАБОТ Методические указания к выполнению практической работы для обучающихся по программам бакалавриата Строительство, Архитектура, Дизайн архитектурной среды, Менеджмент, Технологические машины и оборудование, Ландшафтная архитектура и программам подготовки специалиста...»

«МИНИСТЕРСТВО ОБРАЗОВАНИЯ И НАУКИ РОССИЙСКОЙ ФЕДЕРАЦИИ Федеральное агентство по образованию Государственное образовательное учреждение высшего профессионального образования РОСТОВСКИЙ ГОСУДАРСТВЕННЫЙ СТРОИТЕЛЬНЫЙ УНИВЕРСИТЕТ УТВЕРЖДЕНО на заседании кафедры экономики и управления в строительстве 26 января 2010г. МЕТОДИЧЕСКИЕ УКАЗАНИЯ по выполнению научно-исследовательской работы для студентов, магистрантов и аспирантов экономических специальностей и направлений Ростов-на-Дону, УДК 69.003(07)...»

«База нормативной документации: www.complexdoc.ru Справочные материалы Материалы и технологии века Добромыслов А.Я., Санкова Н.В. Пластмассовые трубы и современные технологии для строительства и ремонта трубопроводов ПРОЕКТИРОВАНИЕ, МОНТАЖ И ЭКСПЛУАТАЦИЯ СИСТЕМ КАнаЛИЗАЦИИ ИЗ ПЛАСТМАССОВЫХ ТРУБ ДЛЯ ЗДАНИЙ И МИКРОРАЙОНОВ РЕКОМЕНДАЦИИ Москва 2004 ПРЕДИСЛОВИЕ Глава 1. ТРУБЫ И ФАСОННЫЕ ЧАСТИ ИЗ ПОЛИМЕРНЫХ МАТЕРИАЛОВ ДЛЯ МОНТАЖА ТРУБОПРОВОДОВ СИСТЕМ КАНАЛИЗАЦИИ ЗДАНИЙ И МИКРОРАЙОНОВ 1.1. Внутренняя...»

«МИНИСТЕРСТВО ОБРАЗОВАНИЯ РОССИЙСКОЙ ФЕДЕРАЦИИ Государственное образовательное учреждение высшего профессионального образования “Оренбургский государственный университет” Кафедра технологии строительных материалов и изделий Т.И. ШЕВЦОВА МАТЕРИАЛОВЕДЕНИЕ МЕТОДИЧЕСКИЕ УКАЗАНИЯ ПО ИЗУЧЕНИЮ ДИСЦИПЛИНЫ “МАТЕРИАЛОВЕДЕНИЕ” Рекомендовано к изданию Редакционно-издательским советом Государственного образовательного учреждения высшего профессионального образования “Оренбургский государственный университет”...»

«Министерство образования и науки Российской Федерации Санкт-Петербургский государственный архитектурно-строительный университет Г. П. КОМИНА, А. О. ПРОШУТИНСКИЙ ГИДРАВЛИЧЕСКИЙ РАСЧЕТ И ПРОЕКТИРОВАНИЕ ГАЗОПРОВОДОВ Учебное пособие Санкт-Петербург 2010 УДК 622.691.4(075.8) Рецензенты: канд. техн. наук, доц. М. А. Кочергин, главный специалист отдела технического надзора Управления капитального строительства ОАО Газпромрегионгаз; А. Г. Матвеев, зам. генерального директора Института...»

«ДЕПАРТАМЕНТ ТРУДА И СОЦИАЛЬНОЙ ПОДДЕРЖКИ НАСЕЛЕНИЯ ЯРОСЛАВСКОЙ ОБЛАСТИ Реализация областной целевой программы Доступная среда. Организация работы органов социальной защиты населения и учреждений социального обслуживания населения Ярославской области по социальной реабилитации инвалидов СБОРНИК ИНФОРМАЦИОННЫХ И МЕТОДИЧЕСКИХ МАТЕРИАЛОВ Ярославль 2011 Реализация областной целевой программы Доступная среда. Организация работы органов социальной защиты населения и учреждений социального обслуживания...»

«МИНИСТЕРСТВО ОБРАЗОВАНИЯ И НАУКИ РОССИЙСКОЙ ФЕДЕРАЦИИ Федеральное государственное бюджетное образовательное учреждение высшего профессионального образования ТЮМЕНСКИЙ ГОСУДАРСТВЕННЫЙ АРХИТЕКТУРНО-СТРОИТЕЛЬНЫЙ УНИВЕРСИТЕТ Кафедра государственного и муниципального управления и права Храмцов А. Б. ГЕОПОЛИТИКА МЕТОДИЧЕСКИЕ УКАЗАНИЯ ПО ПРАКТИЧЕСКИМ ЗАНЯТИЯМ И САМОСТОЯТЕЛЬНОЙ РАБОТЕ СТУДЕНТОВ для направления 081100.62 Государственное и муниципальное управление очной и заочной формы обучения Тюмень,...»

«Федеральное агентство по образованию КАЗАНСКИЙ ГОСУДАРСТВЕННЫЙ АРХИТЕКТУРНО-СТРОИТЕЛЬНЫЙ УНИВЕРСИТЕТ Н.С.Громаков ХИМИЧЕСКАЯ СВЯЗЬ Методические указания по химии для студентов 1 курса дневной, заочной и дистанционной форм обучения Казань 2007 1 УДК 539.19: 541.5(075) ББК 24 К 78 Громаков Н.С. Химическая связь: Методические указания по химии для студентов дневной, заочной и дистанционной форм обучения, Казань: КГАСУ, 2007. -37с. Методические указания содержат основной информационный материал,...»

«МИНИСТЕРСТВО ОБРАЗОВАНИЯ И НАУКИ УКРАИНЫ ХАРЬКОВСКАЯ НАЦИОНАЛЬНАЯ АКАДЕМИЯ ГОРОДСКОГО ХОЗЯЙСТВА В.И. ОСпИщЕВ ОСНОВЫ МАРКЕТИНГА Учебное пособие (для студентов специальности 6.070101 – Транспортные технологии) Харьков Издательство “ФОРТ” 2009 УДК 339.138(075.8) ББК 65.290-2я7 О75 Рецензенты: А.С. Иванилов, д.э.н., профессор, зав. кафедрой экономики Харь­ ковского государственного технического университета строитель­ ства и архитектуры; Г.В. Ковалевский, д.э.н., профессор кафедры маркетинга и ме­...»

«Федеральное агентство по образованию Государственное образовательное учреждение высшего профессионального образования Владимирский государственный университет Кафедра строительных конструкций МЕТОДИЧЕСКИЕ УКАЗАНИЯ К ИЗУЧЕНИЮ РАЗДЕЛА ЖЕЛЕЗОБЕТОННЫЕ КОНСТРУКЦИИ КУРСА КОНСТРУКЦИИ ГОРОДСКИХ ЗДАНИЙ И СООРУЖЕНИЙ Составители: В.В. МИХАЙЛОВ В.И. ВОРОНОВ Владимир 2009 УДК 624.012.3/4 ББК 38.53 М54 Рецензент Кандидат технических наук, доцент зав. кафедрой строительных конструкций Владимирского...»

«МИНИСТЕРСТВО ОБРАЗОВАНИЯ РОССИЙСКОЙ ФЕДЕРАЦИИ Государственное образовательное учреждение высшего профессионального образования – Оренбургский государственный университет А.С. КИЛОВ ОБРАБОТКА МАТЕРИАЛОВ ДАВЛЕНИЕМ В ПРОМЫШЛЕННОСТИ Рекомендовано Ученым советом государственного образовательного учреждения высшего профессионального образования Оренбургский государственный университет в качестве учебного пособия для студентов обучающихся по программам высшего профессионального образования по...»

«Федеральное агентство по образованию Ульяновский государственный технический университет Химия воды Учебное пособие для студентов нехимических специальностей Составители: Л.В. Петрова, Е.Н. Калюкова Ульяновск 2004 УДК 541.1(075.8) ББК 24 Я7 Х 46 Рецензенты: Зав. научно-производственным отделом ООО Трансстройкомплект, канд. техн. наук. И.А. Дорофеев, Зав. кафедрой Химия УлГПУ, доцент, канд. хим. наук И. Т. Гусева Утверждено редакционно-издательским советом университета в качестве учебного...»

«Министерство образования и науки Государственное образовательное учреждение высшего профессионального образования Пермский государственный университет ЕСТЕСТВЕННОНАУЧНЫЙ ИНСТИТУТ Н.Г. Максимович, Е.А. Хайрулина ГЕОХИМИЧЕСКИЕ БАРЬЕРЫ И ОХРАНА ОКРУЖАЮЩЕЙ СРЕДЫ Учебное пособие Пермь 2011 УДК 504.06:550.4 ББК 20.18:26.30 M 18 Максимович, Н.Г. М18 Геохимические барьеры и охрана окружающей среды: учеб. пособие / Н.Г. Максимович, Е.А. Хайрулина; Перм. гос. ун-т. – Пермь, 2011. – 248с.: ил. ISBN...»

«МИНИСТЕРСТВО ОБРАЗОВАНИЯ И НАУКИ УКРАИНЫ ХАРЬКОВСКИЙ НАЦИОНАЛЬНЫЙ УНИВЕРСИТЕТ ГОРОДСКОГО ХОЗЯЙСТВА имени А. Н. БЕКЕТОВА КРОВЕЛЬНЫЕ И ГИДРОИЗОЛЯЦИОННЫЕ РАБОТЫ УЧЕБНОЕ ПОСОБИЕ для студентов высших учебных заведений, которые обучаются по специальности Промышленное и гражданское строительство Под редакцией В. Д. Жван Харьков ХНУГХ 2013 УДК (075) ББК 38.654.3я73-6+38.637я73-6 К83 Авторы: Жван Виктор Денисович – кандидат технических наук, профессор, профессор кафедры Технология...»

« МЕНЕДЖМЕНТ Учебное пособие Для студентов вузов В двух частях Часть 1 Кемерово 2008 2 УДК 65.018 (075) ББК 30.607я7 М 31 Рецензенты: Е.Г. Ягупа, канд. экон. наук, доцент, зав. кафедрой Экономическая теория и экономика предприятий КГСХИ; С.М. Бугрова, канд. экон. наук, доцент кафедры Экономика и организация машиностроительной...»

«Министерство образования и науки Российской Федерации ГОУ ВПО Тамбовский государственный технический университет О.В. УМНОВА, О.В. ЕВДОКИМЦЕВ СТАЛЬНОЙ КАРКАС ЗДАНИЯ ПАВИЛЬОННОГО ТИПА Утверждено Учёным советом университета в качестве учебного пособия Тамбов Издательство ТГТУ 2008 УДК 624.014.2(075) ББК Н549 У545 Р е це н зе н ты: Доктор технических наук, профессор ТГТУ В.И. Леденёв Генеральный директор ООО ФСК Тамбоврегионстрой В.И. Скрылев Умнова, О.В. У545 Стальной каркас здания павильонного...»

«Министерство образования и науки Российской Федерации ФГБОУ ВПО Казанский государственный архитектурно-строительный университет Кафедра Экономики и предпринимательства в строительстве Методические указания Организация презентаций. Использование мультимедиа-ресурсов для студентов следующих специальностей и направлений подготовки: 070603 Искусство интерьера, 190702 Организация и безопасность движения, 190205 Подъемно-транспортные, строительные, дорожные машины и оборудование, 270100.62...»

«МИНИСТЕРСТВО ОБРАЗОВАНИЯ И НАУКИ РОССИЙСКОЙ ФЕДЕРАЦИИ ФЕДЕРАЛЬНОЕ АГЕНТСТВО ПО ОБРАЗОВАНИЮ САНКТ-ПЕТЕРБУРГСКИЙ НАЦИОНАЛЬНЫЙ ИССЛЕДОВАТЕЛЬСКИЙ УНИВЕРСИТЕТ ИНФОРМАЦИОННЫХ ТЕХНОЛОГИЙ, МЕХАНИКИ И ОПТИКИ А.А. Воробьева УЧЕБНОЕ ПОСОБИЕ ПО КУРСУ ГЕОИНФОРМАЦИОННЫЕ СИСТЕМЫ ТЕРРИТОРИАЛЬНОГО УПРАВЛЕНИЯ Санкт-Петербург 2012 1 Учебное пособие посвящено геопространственному моделированию объектов с помощью ГИС и использование сопровождаемой их семантической информации. Кроме того вопросам сбора и подготовки...»

«МИНИСТЕРСТВО ОБРАЗОВАНИЯ И НАУКИ РОССИЙСКОЙ ФЕДЕРАЦИИ КАЗАНСКИЙ ГОСУДАРСТВЕННЫЙ АРХИТЕКТУРНО-СТРОИТЕЛЬНЫЙ УНИВЕРСИТЕТ Кафедра производственной безопасности и права БЕЗОПАСНОСТЬ ЖИЗНЕДЕЯТЕЛЬНОСТИ Методические указания для выполнения контрольной работы для студентов-заочников направления Строительство специальности 270106.65 и профилю 270804.62 Производство и применение строительных материалов, изделий и конструкций Казань 2013 УДК 69.05: 658.382 ББК К 66 К 66 Безопасность жизнедеятельности:...»

А.В. Игнатьев, Н.А. Михайлова, Т.В. Ерещенко МЕТОД КОНЕЧНЫХ ЭЛЕМЕНТОВ И ЕГО РЕАЛИЗАЦИЯ В СРЕДЕ MATHCAD Лабораторный практикум по дисциплине «Аналитические и численные методы решения уравнений математической физики» Волгоград 2010 Министерство образования и науки Российской Федерации Волгоградский государственный архитектурно-строительный университет Кафедра прикладной математики и вычислительной техники А.В. Игнатьев, Н.А. Михайлова, Т.В. Ерещенко МЕТОД КОНЕЧНЫХ ЭЛЕМЕНТОВ И ЕГО РЕАЛИЗАЦИЯ В СРЕДЕ MATHCAD Лабораторный практикум по дисциплине «Аналитические и численные методы решения уравнений математической физики» Волгоград 2010 УДК 624.04:004.92(076.5) ББК 38.112я73+32.973-018.2я73 И 266 Р е ц е н з е н т ы: кандидат технических наук М.М. Степанов, профессор кафедры прикладной математики и вычислительной техники ВолгГАСУ; доктор технических наук Н.Г. Бандурин, профессор кафедры строительной механики ВолгГАСУ Утверждено редакционно-издательским советом университета в качестве учебно-практического пособия Игнатьев А.В. И 266 Метод конечных элементов и его реализация в среде Mathcad: лабораторный практикум по дисциплине «Аналитические и численные методы решения уравнений математической физики» / А.В. Игнатьев, Н.А. Михайлова, Т.В. Ерещенко; Волгогр. гос. архит.строит. ун-т. Волгоград: ВолгГАСУ, 2010. 31 с. ISBN 978-5-98276-372-3 Содержатся краткие теоретические сведения, необходимые для выполнения лабораторной работы по дисциплине «Аналитические и численные методы решения уравнений математической физики», приведены варианты индивидуальных заданий, примеры выполнения заданий лабораторной работы, а также сформулированы контрольные вопросы по изучаемой теме. Для магистров специальности АД, СМ и ВиВ дневной формы обучения. УДК 624.04:004.92(076.5) ББК 38.112я73+32.973-018.2я73 ISBN 978-5-98276-372-3 Государственное образовательное учреждение высшего профессионального образования «Волгоградский государственный архитектурно-строительный университет», 2010 2 Лабораторная работа. МЕТОД КОНЕЧНЫХ ЭЛЕМЕНТОВ И ЕГО РЕАЛИЗАЦИЯ В СИСТЕМЕ MATHCAD Цель работы: изучить метод конечных элементов и получить навыки его реализации в системе Mathcad. Основные понятия и концепция МКЭ Основная идея метода Основная идея метода состоит в представлении рассчитываемой конструкции в виде совокупности элементов простой формы, соединенных между собой в отдельных точках. По сути дела, сплошная среда с бесконечным числом степеней свободы заменяется набором подобластей, имеющих конечное число степеней свободы. При таком подходе искомые непрерывные величины (перемещения, напряжения, деформации и т. д.) внутри каждого конечного элемента (КЭ) выражаются с помощью аппроксимирующих функций через узловые значения этих величин. Распределенные внешние нагрузки заменяются эквивалентными узловыми силами. В математическом плане задача состоит в приведении дифференциальных уравнений или энергетического функционала, описывающих рассматриваемую конструкцию, к системе алгебраических уравнений, решение которой дает значения искомых узловых неизвестных. Метод конечных элементов имеет очень широкое распространение в практике расчетов на прочность, устойчивость и колебания строительных, машиностроительных, авиационных конструкций. С помощью МКЭ можно успешно выполнить анализ широкого класса стержневых систем (фермы, рамы и т. д.), тонкостенных пространственных конструкций (плиты перекрытий, оболочки покрытий и т. д.), массивных трехмерных тел, а также комбинированных систем, состоящих из одномерных, двумерных и трехмерных элементов. МКЭ отличает широкая область применимости, инвариантность по отношению к геометрии конструкции и физическим характеристикам материалов, относительная простота учета взаимодействия конструкций с окружающей средой (механические, температурные, коррозионные воздействия, граничные условия и т. д.), высокая степень приспособляемости к автоматизации всех этапов расчета. Метод имеет простую физическую интерпретацию и тесно связан с методом перемещений, который широко используется в строительной механике. 3 На базе конечно-элементного подхода разработано большое количество мощных программных комплексов. Среди них можно отметить такие, как ABAQUS, ADINA, ANSYS, COSMOS/М, MSC/NASTRAN, ЛИРА, SCAD, STARK, СТАДИО. Большинство из них имеет обширную библиотеку конечных элементов и дает возможность выполнять расчеты на прочность, устойчивость и колебания, учитывать физическую и геометрическую нелинейности, ортотропию материала, температурные нагрузки и т. д. Представленный выше перечень программных продуктов, реализующих МКЭ, является далеко не полным и лишь отражает ситуацию в данной области в настоящий момент. Несомненно, МКЭ имеет существенные преимущества по сравнению с другими подходами и в значительной степени универсален. Вместе с тем его следует рассматривать как одну из многочисленных ступеней развития средств численного исследования при проектировании. Общая схема алгоритма МКЭ Метод конечных элементов предусматривает следующие основные этапы: 1. Идеализация физической системы. Под идеализацией понимают процесс перехода от исходной физической системы к математической модели. Этот процесс является наиболее важным шагом при решении технической или инженерной задачи. Ключевым пунктом в этом процессе является понятие модели, которую можно определить как символическое устройство, построенное для моделирования и предсказания поведения системы. Математическое моделирование, или идеализация, есть процесс, с помощью которого инженер переходит от реальной физической системы к математической модели системы. Данный процесс называется идеализацией, поскольку математическую модель необходимо абстрагировать от физической реальности. В качестве примера реальной физической системы рассмотрим инженерную конструкцию в виде плоской пластины, нагруженную поперечными силами. Математические модели данной системы, которые инженер может использовать для анализа напряжений в пластине, могут быть следующими: 1) модель очень тонкой пластины, основанная на теории изгиба мембран; 4 2) модель тонкой пластины, основанная на классической теории Кирхгоффа; 3) модель достаточно толстой пластины, основанная, например, на теории Миндлина-Рейсснера; 4) модель очень толстой пластины, основанная на трехмерной теории упругости. Очевидно, инженер должен обладать достаточными теоретическими знаниями, чтобы правильно выбрать соответствующую математическую модель системы (конструкции), которую ему необходимо исследовать. 2. Дискретизация рассматриваемой области. Рассчитываемая конструкция разбивается воображаемыми точками, линиями или поверхностями на элементы конечных размеров (конечные элементы). Предполагается, что элементы связаны между собой в узловых точках, расположенных на их границах. В некоторых задачах строительной механики в качестве искомых неизвестных помимо узловых перемещений принимаются также их частные производные. 3. Построение интерполирующих функций. Выбирается система функций (чаще всего кусочно-полиномиальная), однозначно определяющая перемещения внутри каждого конечного элемента через перемещения узловых точек. Интерполирующие функции подбираются таким образом, чтобы обеспечить непрерывность искомых величин (перемещений и их производных) вдоль границ элемента. 4. Вывод основных геометрических и физических соотношений. На основе выбранной системы интерполирующих функций выводятся зависимости между деформациями и перемещениями (геометрические соотношения), а также между напряжениями и деформациями (физические соотношения). 5. Построение матрицы жесткости конечного элемента. С помощью принципа Лагранжа на основе полученных геометрических и физических соотношений строится матрица жесткости конечного элемента. 6. Получение системы уравнений метода конечных элементов. Каждая матрица жесткости отдельного конечного элемента включается в глобальную матрицу жесткости в цикле по элементам. Таким образом формируется система алгебраических уравнений всей конструкции (уравнения равновесия), которая имеет вид Kz = P , 5 где K - матрица жесткости системы (ансамбля) конечных элементов; z - вектор неизвестных узловых перемещений; Р - вектор узловых нагрузок. В матрице жесткости K, записанной выше системы уравнений, необходимо учесть граничные условия, так как в противном случае эта матрица будет вырожденной. 7. Решение системы алгебраических уравнений. Для решения системы линейных алгебраических уравнений (СЛАУ) используются как точные, так и (при высоком порядке системы) итерационные методы. Построенные на их основе эффективные численные процедуры учитывают симметрию и ленточную структуру матрицы жесткости системы. 8. Определение деформаций и напряжений. Деформации, напряжения и усилия в конструкции определяются с помощью найденных узловых перемещений на основе геометрических и физических соотношений. Рассмотрим некоторые из этих этапов более подробно. Дискретизация рассматриваемой области Разбиение конструкции на конечные элементы - это весьма важный шаг в процедуре решения задачи по МКЭ, поскольку от него во многом зависит точность получаемого решения. Успех на этом этапе обеспечивают, в первую очередь, имеющиеся инженерные навыки. Неудачно выполненное разбиение области на конечные элементы может привести к ошибочным результатам. При назначении сетки КЭ возникает задача оптимального разбиения конструкции на подобласти. При этом следует иметь в виду, что размеры элементов должны быть достаточно малыми для того, чтобы обеспечить приемлемую точность решения, с другой стороны, использование густой сетки приводит к большим системам алгебраических уравнений, решение которых связано со значительным объемом вычислительной работы. В процессе разбиения области на конечные элементы необходимо учитывать некоторые общие представления об окончательных результатах расчета для того, чтобы уменьшить размеры конечных элементов в зонах концентрации напряжений, где искомые величины быстро меняются, и увеличить размеры КЭ там, где искомые величины меняются медленно. Важным моментом в процессе решения задачи по МКЭ является нумерация узлов сетки, поскольку от этого зависит ширина ленты 6 матрицы разрешающих уравнений, соответственно время счета и объем используемой памяти ЭВМ. В настоящее время разработаны сервисные программы автоматизированной разбивки конструкции на конечные элементы и рациональной нумерации узлов. Построение интерполирующих функций МКЭ основан на аппроксимации непрерывной функции, определенной на всей области, дискретной моделью с помощью кусочно-непрерывных функций, определенных на подобластях (конечных элементах). Запишем перемещения, являющиеся функциями координат произвольной точки конечного элемента, через компоненты вектора узловых перемещений с помощью интерполирующей функции (функции формы или базисной функции): u = Nz , (1) где N = [ N1 N 2 … N s ] - матрица функции формы; z = { z1 z2 … zs } - вектор узловых перемещений конечного элемента (КЭ); s - количество степеней свободы КЭ. Функции (1) должны удовлетворять критериям полноты и совместности. Рассмотрим их. 1. Критерий полноты. Интерполирующая функция должна обеспечивать постоянные значения рассматриваемых величин при уменьшении размеров элемента. Для выполнения этого условия интерполирующая функция должна представлять собой полный полином как минимум степени р, где р - наивысший порядок производной, входящей в функционал. T Условие полноты удовлетворяется в том случае, когда s ∑ Ni = 1 . i =1 2. Критерий совместности. Интерполирующая функция должна быть непрерывна вместе со своими производными до (n – 1)-гo порядка включительно (где n - максимальный порядок производной в подынтегральном выражении функционала энергии) на границе между элементами. Критерии полноты и совместности представляют собой достаточные условия сходимости метода конечных элементов. При их выполнении с уменьшением размеров конечного элемента приближенные 7 решения МКЭ монотонно сходятся к точному решению. Сказанное отнюдь не означает, что нарушение данных критериев приводит к невозможности получения достоверного решения. Существуют несовместные и даже неполные элементы, которые дают высокую точность и быструю сходимость. Вывод основных геометрических и физических соотношений В общем виде зависимость между деформациями и перемещениями (геометрические соотношения) записывается так: ε = Bz , (2) где ε - вектор деформаций; z - вектор узловых перемещений; В - матрица, связывающая вектор узловых перемещений с вектором, содержащим компоненты тензора деформации. Таким образом, полагается, что зависимость (2) между деформациями и узловыми перемещениями носит линейный характер. Линейной зависимости отвечают такие условия работы конструкции, когда деформации и углы поворота малы по сравнению с единицей, а квадраты углов поворота малы по сравнению с соответствующими компонентами деформации. Физические соотношения, определяющие зависимость между напряжениями и деформациями, имеют вид σ = Dε, (3) где σ - вектор, содержащий компоненты тензора напряжений; D - матрица упругости. Уравнения состояния (3) представляют собой обобщенный закон Гука, устанавливающий прямо пропорциональную зависимость между напряжениями и деформациями, справедливый для определенного класса материалов на определенном участке графика σ − ε . Построение матрицы жесткости конечного элемента Решение задач расчета конструкций базируется на двух основных подходах. В первом случае решаются дифференциальные уравнения с заданными граничными условиями. Во втором случае записывается условие стационарности интегральной величины, связанной с работой напряжений и внешней приложенной нагрузки и представляющей собой полную потенциальную энергию системы. Для расчета конструкций в рамках МКЭ используется второй подход. Как известно, полная потенциальная энергия упругой системы определяется по формуле 8 Π (z) = W (z) − A(z) , где W - потенциальная энергия деформации; А - потенциал внешних сил. Потенциальная энергия деформации упругой системы определяется соотношением W= 1 T ∫ ε σ dV , 2V где V - объем, занимаемый телом, а потенциал внешних распределенных нагрузок определяется по формуле A = ∫ uT p dS , S где р - вектор внешних распределенных нагрузок; S - площадь, на которой приложена нагрузка. При этом деформации и напряжения, входящие в формулу для потенциальной энергии, выражаются через узловые перемещения. Получение уравнений МКЭ в перемещениях основано на одном из фундаментальных энергетических принципов механики - принципе Лагранжа, согласно которому для системы, находящейся в состоянии равновесия, полная потенциальная энергия принимает стационарное значение. Это условие записывается в виде ∂Π = 0. ∂z Будем считать, что значение полной потенциальной энергии для всей области V равно сумме энергий отдельных конечных элементов: m () m (() ()) Π (z) = ∑ Π z = ∑ W i z i − Ai z i , i =1 i i i =1 (4) где m - количество конечных элементов. Тогда ∂Π m ⎛ ∂W i (z) ∂Ai (z) ⎞ = ∑⎜ − ⎟ = 0. ∂z i =1 ⎝ ∂z ∂z ⎠ (5) Рассмотрим отдельный конечный элемент, индекс i при этом опускаем: 9 1 T 1 T (Bz)T DBz dV − ∫ (Nz)T p dS = ε σ dV − u p dS = ∫ ∫ ∫ 2V 2V S S 1 1 = zT (∫ BT DBz dV) z − zT ∫ N T p dS = zT Kz − zT P, (6) 2 2 S Π(z) = где K = ∫ BT DB dV - матрица жесткости конечного элемента, а (7) V T P = ∫ N p dS - вектор узловых нагрузок. (8) S Получение системы уравнений метода конечных элементов Для выполнения операции суммирования необходимо преобразовать векторы узловых перемещений z(i) и узловой нагрузки P(i) для отдельного конечного элемента в соответствующие векторы z и Р для всей системы, что может быть сделано с помощью некоторой булевой матрицы H(i), содержащей в качестве элементов только нули и единицы: z (i) = H (i) z; P (i) = H (i) P. (9) Подстановка формул (9) в выражение для полной потенциальной энергии конечного элемента (6) дает: () () T 1 (i) T (i) (i) H z K H z − H (i) z H (i) P = 2 T T 1 = zT H (i) K (i) H (i) z − zT H (i) H (i) P. 2 Тогда дифференцирование по z, в соответствии с формулой (5), приводит к системе уравнений: Π (i) = m ∑ i =1 (T T) H (i) K (i) H (i) z − H (i) H (i) P = 0 , (10) где использовано правило дифференцирования матричных соотношений ∂ T z Kz = 2 Kz . ∂z Система (10), которую можно записать в виде () Kz = P , (11) 10 представляет собой систему линейных алгебраических уравнений метода конечных элементов, являющихся уравнениями равновесия в перемещениях. Как правило, решение системы (11) выполняется методом Гаусса или итерационными методами. (i)T Матрица жесткости отдельного элемента H K (i) H (i) , фигурирующая в формуле (10), представляет собой расширенную матрицу, размерность которой равна размерности глобальной матрицы. Поэтому использование процедуры суммирования в формуле (10) при численной реализации МКЭ неэффективно. В практических расчетах выполняется прямое построение глобальной матрицы жесткости. В этом случае строится матрица K для отдельного конечного элемента по формуле (7), имеющая размерность S×S. Затем строкам и столбцам этой матрицы приписываются номера глобальных степеней свободы, что позволяет определить местоположение коэффициентов матрицы жесткости КЭ в глобальной матрице жесткости. После этого в предварительно обнуленную глобальную матрицу заносятся коэффициенты матрицы жесткости КЭ на то место, которое определено их адресом. Пусть, например, система состоит из двух КЭ, содержащих по два узла, в каждом из которых имеется одно неизвестное (одна степень свободы). Общее количество узлов системы - 3, размерность глобальной матрицы 3×3, элементы соединены между собой во 2-м узле. Матрицы жесткости 1-го и 2-го КЭ, с соответствующей нумерацией коэффициентов, и глобальная матрица имеют вид K (1) = 1 ⎡ k11 ⎢ 1 ⎣⎢ k21 1 ⎤ k12 ⎥; 1 k22 ⎦⎥ K (2) = 2 ⎡ k22 ⎢ 2 ⎢⎣ k32 1 1 ⎡ k11 k12 2 ⎤ ⎢ 1 k23 1 2 = ; K ⎥ ⎢ k21 k22 + k22 2 k33 ⎥⎦ ⎢ 2 k32 ⎢⎣ 0 0 ⎤ ⎥ 2 k23 ⎥. 2 ⎥ k33 ⎥⎦ В матрице K системы уравнений (11) необходимо учесть граничные условия, в противном случае она будет вырожденной, т. е. ее определитель будет равен нулю. Учет граничных условий может быть осуществлен тремя различными способами 1. Из матрицы K удаляются k-я строка и k-й столбец, соответствующие перемещению zk = 0 . После этого строки и столбцы матрицы 11 перенумеровываются. Соответственно уменьшается размерность вектора узловых перемещений. 2. Уравнение zk = 0 , соответствующее граничному условию, формируется в составе матрицы K. Для получения zk = 0 в матрице K k-я строка и k-й столбец, а также соответствующий элемент в векторе внешних нагрузок P заполняются нулями. На место диагонального элемента rrr в матрице K ставится единица. В результате порядок матрицы не изменяется, а заданные перемещения получают нулевые значения. 3. Для получения zk = 0 диагональный элемент rrr умножается на большое число. Порядок матрицы при этом не меняется. Определение деформаций и напряжений В результате решения системы уравнений (11) определяется вектор узловых перемещений всей конструкции. На основе найденных значений узловых перемещений по формуле (2) определяется вектор деформаций КЭ, а по формуле (3) - вектор напряжений. Двумерный симплекс-элемент Классификация конечных элементов может быть проведена в соответствии с порядком многочленов - функций этих элементов. При этом рассматриваются три группы элементов: 1) симплекс-элементы; 2) комплекс-элементы; 3) мультиплекс-элементы. Симплекс-элементам соответствуют многочлены первой степени. Комплекс-элементам - многочлены более высокого порядка. В симплекс-элементе число узлов равно размерности пространства +1. В комплекс-элементе число узлов больше этой величины. Для мультиплекс-элементов также используются многочлены высокого порядка, но границы элементов при этом должны быть параллельны координатным осям. Рассмотрим формирование матрицы жесткости для двумерного симплекс-элемента. Двумерный симплекс-элемент представляет собой треугольник с узлами, расположенными в его вершинах (рис. 1). 12 Рис. 1. Двумерный симплекс-элемент Узлы КЭ нумеруются против часовой стрелки, начиная с некоторого произвольно выбранного i-го узла. Координаты i-гo, j-го и k-гo узлов по оси x обозначены через xi , x j , xk , по оси y - через yi , y j , yk . Каждый узел имеет две степени свободы - перемещение u вдоль оси x и перемещение v вдоль оси y. Интерполирующие функции, определяющие перемещение произвольной точки КЭ вдоль осей х и у, принимаются в виде u (x, y) = α 0 + α1 x + α 2 y; (12) v(x, y) = α3 + α 4 x + α5 y. Коэффициенты α 0 ,…, α5 определяются с помощью граничных условий: u = ui , v = vi при x = xi и y = yi ; u = u j , v = v j при x = x j и y = y j ; u = uk , v = vk при x = xk и y = yk . Определим коэффициенты α 0 , α1 , α 2 . Для этого подставим граничные условия для функции и в первое выражение (12), что приведет к системе уравнений: ui = α 0 + α1 xi + α 2 yi ; u j = α 0 + α1 x j + α 2 y j ; uk = α 0 + α1 xk + α 2 yk . 13 ⎡1 xi ⎢ Или ⎢1 x j ⎢1 x k ⎣ yi ⎤ ⎧α 0 ⎫ ⎧ ui ⎫ ⎥⎪ ⎪ ⎪ ⎪ y j ⎥ ⎨ α1 ⎬ = ⎨u j ⎬ . ⎪ ⎪ ⎪ ⎪ yk ⎥⎦ ⎩α 2 ⎭ ⎩uk ⎭ (13) Определитель системы (13) равен удвоенной площади F треугольного элемента: 1 xi 1 xj yi y j = 2F . 1 xk yk (14) Тогда по правилу Крамера α0 = ui uj xi xj yi yj uk xk yk 1 xi 1 xj yi yj 1 xk yk = ui uj xi xj yi yj uk xk yk 2F или 1 ⎡ ui x j yk − x k y j − u j (xi yk − x k yi) − uk xi y j − x j yi ⎤ . ⎣ ⎦ 2F Аналогично определяются коэффициенты α1 и α 2 . После подстановки выражений для α0 , α1 , α 2 в первую формулу (12) имеем 1 ⎡ u(x, y) = (ai + bi x + c i y) ui + a j + bj x + c j y u j + 2F ⎣ + (ak + bk x + c k y) uk ⎤⎦ , (15) α0 = () (()) где ai = x j yk − x k y j ; bi = y j − yk ; ci = xk − x j . (16) Остальные коэффициенты в (15) получаются по формулам (16) циклической перестановкой индексов (индекс i заменяется на индекс j, индекс j - на индекс k, индекс k - на индекс i). Аналогично: v(x, y) = 1 ⎡ (ai + bi x + c i y) vi + a j + bj x + c j y v j + 2F ⎣ + (ak + bk x + c k y) vk ⎤⎦ . () 14 (17) Тогда в матричном виде ⎡ Ni ⎧u ⎫ u = ⎨ ⎬ = Nz = ⎢ ⎩v ⎭ ⎢⎣ 0 0 Nj 0 Nk Ni 0 Nj 0 ⎧ ui ⎫ ⎪v ⎪ ⎪ i⎪ 0 ⎤ ⎪⎪u j ⎪⎪ ⎥ ⎨ ⎬, N k ⎥⎦ ⎪ v j ⎪ ⎪u ⎪ ⎪ k⎪ ⎩⎪ vk ⎭⎪ 1 1 a j + bj x + c j y ; (ai + bi x + c i y) ; N j = 2F 2F (18) 1 Nk = (ak + bk x + c k y). 2F Геометрические соотношения, связывающие деформации и перемещения в рамках плоской задачи теории упругости, записываются с использованием формул (15), (17) следующим образом: (где Ni =) ∂u 1 ∂v 1 bi ui + b j u j + bk uk ; ε y = ci vi + c j v j + ck vk ; = = ∂x 2 F ∂y 2 F ∂u ∂v 1 ci ui + c j u j + ck uk + bi vi + b j v j + bk vk γ xy = + = ∂y ∂x 2 F или в матричном виде: εx = () ((⎧ε ⎫ ⎡bi x ⎪ ⎪ 1 ⎢ ε = ⎨ εy ⎬ = ⎢0 ⎪ ⎪ 2F ⎢ ⎩ γ xy ⎭ ⎣ci ⎡bi 1 ⎢ где B = ⎢0 2F ⎢ ⎣ci)) 0 bj 0 bk ci 0 cj 0 bi cj bj ck 0 bj 0 bk ci 0 cj 0 bi cj bj ck ⎧ ui ⎫ ⎪v ⎪ 0 ⎤⎪ i ⎪ ⎥ ⎪⎪u j ⎪⎪ ck ⎥ ⎨ ⎬ = Bz , ⎥ ⎪v j ⎪ bk ⎦ ⎪ ⎪ u ⎪ k⎪ ⎩⎪ vk ⎭⎪ 0⎤ ⎥ ck ⎥ - матрица градиентов. ⎥ bk ⎦ 15 (19) Величина удвоенной площади конечного элемента 2F в выражении (19) подсчитывается по формуле (14). Физические соотношения, определяющие зависимость между напряжениями и деформациями в плоской задаче теории упругости, записываются в виде σ = Dε , ⎡ ⎤ ⎢ 1 ν1 0 ⎥ ⎥ E1 ⎢ где D = 1 0 ν 1 ⎥ - матрица упругости. 2 ⎢ 1 − ν1 ⎢ 1 − ν1 ⎥ ⎢0 0 ⎥ 2 ⎦ ⎣ (20) В случае плоской деформации (ε z = 0) в формуле (20) следует принять E1 = E ν ; ν = , 1 1− ν 1 − ν12 (21) а для обобщенного плоского напряженного состояния (σ z = 0) E1 = E ; ν1 = ν. (22) Формулы (21) и (22) соответствуют изотропному материалу с модулем упругости Е и коэффициентом Пуассона ν . Не представляет большого труда построить матрицы упругости для ортотропного материала, когда жесткостные характеристики различны в двух взаимно перпендикулярных направлениях. Поскольку матрицы В и D содержат только константы, объемный интеграл, определяющий матрицу жесткости элемента в формуле (7), легко вычисляется: K = ∫ BT DB dV = BT DB ∫ dV V (23) V или K = BT DB tF . (24) В формуле (24) t - толщина элемента; F - площадь элемента. Обычно матрица жесткости (24) определяется численно. Для этого вначале находятся численные значения коэффициентов мат16 риц В и D, а затем выполняется перемножение в соответствии с выражением (23) или (24). ПРИМЕР ВЫПОЛНЕНИЯ ЛАБОРАТОРНОЙ РАБОТЫ Перед выполнением лабораторной работы еще раз напомним основные этапы метода МКЭ: 1. Упругое тело разбивается на элементы. Объемное тело на тетраэдры или параллелепипеды. Плоское тело - на треугольники и прямоугольники. 2. Для каждого элемента составляется матрица жесткости с использованием функции формы. Функция формы представляет собой способ аппроксимации неизвестной функции перемещений. 3. Матрицы жесткости элементов объединяются в единую матрицу жесткости для всего тела. 4. Решая систему уравнений, находят узловые перемещения. 5. С помощью уравнений теории упругости определяются деформации и напряжения в узловых точках тела. В приведенном примере решается плоская задача теории упругости. Кольцо, нагруженное двумя силами (рис. 2, в), имеет две оси симметрии, поэтому для повышения точности расчета рассмотрим одну четвертую часть кольца (рис. 3). На осях симметрии должны выполняться граничные условия равенства нулю перемещений, перпендикулярных осям симметрии. Рассматриваемую четверть кольца разбиваем на треугольные конечные элементы (см. рис. 2, б). Треугольный элемент имеет 6 степеней свободы (независимых узловых перемещений). Нумерация узловых перемещений в элементе начинается в нижнем левом узле треугольника и продолжается против часовой стрелки. Горизонтальные перемещения - нечетные, вертикальные - четные. Нумерация узлов всего тела и конечных элементов - по столбцам сверху вниз, слева направо. Размеры элементов могут быть разными (чем меньше элемент, тем выше точность расчетов). В нашем примере всего 66 узлов и 100 конечных элементов. Положение рассчитанных узлов показано на рис. 3, б. Расчет координат узлов приведен на рис. 4. Ввод координат узлов - это кропотливый труд, и при большом количестве узлов лучше автоматизировать эту работу. 17 а б в Рис. 2. Схема нагружения и треугольный конечный элемент а б Рис. 3. Расчетная схема и координаты узлов На рис. 4 приводится расчет полярных координат узлов и их преобразование в прямоугольные (декартовы) координаты. Здесь r1 и r2 - наружный и внутренний радиусы кольца; t - толщина кольца; φ1 и φ2 - начальное и конечное значения угловой координаты; X0 и Y0 - декартовы координаты полюса (начала полярных координат); nr и nφ - число узлов в столбце (вдоль радиуса) и в ряду (по углу охвата рассматриваемой части тела). Результаты расчетов координат узлов приведены на графике (рис. 5). 18 Рис. 4. Расчет координат узлов элементов Рис. 5. Сетка узлов 19 Кропотливой и трудоемкой является также задача составления матрицы индексов. На рис. 6 приведена используемая в примере программа составления матрицы индексов. Здесь же приведен автоматический расчет граничных условий и, в зависимости от числа узлов, номеров перемещений, в которых перемещение на осях симметрии равно нулю. Автоматизация расчетов координат узлов, матрицы индексов и граничных условий позволяет для данной схемы менять количество узлов (рис. 7). Рис. 6. Программа расчета матрицы индексов 20 Рис. 7. Матрица индексов, координаты узлов, номера и заданные перемещения Координаты узлов, матрицу индексов и граничные условия можно записать в отдельные файлы, как это показано на рис. 8. В дальнейшем эти файлы можно считать с помощью функции READPRN и использовать в другом документе. 21 Рис. 8. Запись координат, матрицы индексов и граничных условий во внешние файлы Данный расчет проведен с учетом размерностей (рис. 4). Учет размерностей вносит дополнительные трудности в достаточно сложный расчет, особенно при вводе матриц от размерных величин (рис. 9). Рис. 9. Вектор сил и матрица индексов перемещений 22 У каждого треугольного элемента 3 узла и 6 узловых перемещений. Матрица индексов перемещений (матрица связи глобальных номеров узловых перемещений тела с локальными номерами узловых перемещений элементов) получена путем удвоения матрицы MIU. Матрица жесткости элемента рассчитывается по формуле K = ∫ BT DB dV = BT DB ∫ dV . V (23) V Здесь B = ∂T N , где D - матрица внутренней жесткости, содержащая упругие постоянные материала E , ν ; ∂ - матричный дифференциальный оператор, означающий определенную последовательность присвоения знака дифференцирования; N - матрица функций формы. Для треугольного элемента функция формы - уравнение плоскости, определяемое выражениями (18). Как было показано выше (см. выражение (19)), матрица B = ∂T N содержит константы, которые зависят только от координат узлов. Приведем расчет коэффициентов для формирования матрицы жесткости элемента (рис. 10) и расчет площади элементов (рис. 11). Рис. 10. Расчет коэффициентов для формирования матрицы жесткости элемента 23 Рис. 11. Расчет площади элементов Матрица внутренней жесткости D приведена ниже (рис. 12) и записана в виде условного оператора: разные матрицы для плоского напряженного NDS = 0 и плоского деформированного состояния NDS = 1. Рис. 12. Формирование матрицы внутренней жесткости Для треугольного элемента интеграл по объему равен произведению подынтегрального выражения на объем. Формула (23) для расчета матрицы жесткости элемента приведена вверху. Матрица жесткости системы формируется с помощью матрицы индексов. 24 Учет граничных условий сопровождается перестройкой матрицы жесткости системы и вектора сил (рис. 13). Рис. 13. Формирование матрицы жесткости системы и учет граничных условий Узловые перемещения определяются путем обращения матрицы жесткости. Рис. 14. Определение перемещений узлов, деформаций и напряжений в центре элемента Согласно уравнениям теории упругости ε = ∂T u , где u - вектор перемещений. Согласно уравнениям связи узловых перемещений ∆ и перемещений u произвольной точки u = N ∆ . 25 () Отсюда деформация элемента ε = ∂T N ∆ . Из физических уравнений теории упругости (закон Гука) напряжения σ = Dε . Сложность расчета состоит в аккуратном использовании индексов элементов, узлов, столбцов, строк, присвоении индексам значений, взятых из матрицы индексов. Для треугольного элемента функция формы линейна, поэтому производные от функции формы, деформации и напряжения, найденные на рис. 14, постоянны во всей площади элемента. Напряжения в узлах тела определяются как среднее арифметическое напряжений или деформаций во всех элементах, сходящихся в узле. Расчет напряжений и деформаций в узлах тела приведен на рис. 15. Рис. 15. Определение перемещений узлов тела, напряжений и деформаций в центре каждого элемента 26 На рис. 15 показано определение 4-го значения деформации и напряжения в каждом элементе, не учитываемых в матрице внутренней жесткости D. Работая с документом Mathcad, если мы опустим выражение NDS: = 1 ниже выражения NDS: = 0, то мы можем увидеть результаты расчета уже не при плоском напряженном состоянии, а при плоской деформации. Результаты расчета приведены на рис. 16. Рис. 16. Результаты расчета 27 Задание к лабораторной работе Решить методом конечных элементов (МКЭ) плоскую задачу теории упругости. Кольцо, нагруженное двумя силами (рис. 3, б), имеет две оси симметрии, поэтому для повышения точности расчета необходимо рассмотреть одну четвертую часть кольца. При помощи подходящей сетки она разделяется на систему треугольных конечных элементов. Число узлов вдоль радиуса nr и число узлов по углу охвата рассматриваемой части тела nφ выбираются из таблицы индивидуальных заданий в соответствии с вариантом. Определить деформации и напряжения при плоском напряженном состоянии и при плоском деформированном состоянии. Варианты индивидуальных заданий № Число узлов вари- вдоль радиуанта са, nr 1 8 2 9 3 10 4 11 5 12 6 13 7 8 8 9 9 10 10 11 11 12 12 13 13 8 14 9 15 10 Число узлов по углу охвата, nφ 4 5 6 7 8 9 5 6 7 8 9 4 6 7 8 № варианта 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 Число узлов вдоль радиуса, nr 11 12 13 9 10 11 12 13 10 11 12 13 11 12 13 Число узлов по углу охвата, nφ 9 4 5 6 4 8 9 7 5 6 7 6 4 5 6 Содержание отчета В рабочем каталоге студента должны быть созданы два файла, содержащие отлаженные документы в системе Mathcad, соответствующие расчетам при плоском напряженном состоянии и при плоском деформированном состоянии. 28 Отчет по лабораторной работе должен содержать: 1) название лабораторной работы; 2) цель лабораторной работы; 3) задание; 4) скопированные с экрана монитора отлаженные документы выполнения задания. Контрольные вопросы 1. Для чего предназначен метод конечных элементов? 2. Привести общую схему расчета по методу конечных элементов. 3. Что такое идеализированная конструкция? 4. Из какой системы линейных алгебраических уравнений определяются перемещения узлов? 5. Как определяется вид аппроксимирующего полинома? 6. Какие функции называются базисными функциями? 7. Как определяются функции формы? 8. Как определяются функции перемещений? 9. Как определяются функции деформаций? 10. Как определяются напряжения? 11. Какой принцип используется для определения обобщенных сил? 12. Как определяется матрица жесткости элемента? 13. Как составляется матрица жесткости системы? 14. От учета каких воздействий получаются свободные члены системы канонических уравнений? 15. Чем могут быть вызваны предварительные деформации? 16. Чем могут быть вызваны предварительные напряжения? 17. Как определяются реактивные усилия вследствие отдельных воздействий? 18. Как в системе Mathcad производится расчет координат узлов сетки? 29 Библиографический список 1. Дарков В.А. Строительная механика: учеб. для строит. спец. вузов / В.А. Дарков, Н.Н. Шапошников. Изд. 8-е, перераб. и доп. М. : Высш. шк., 1986. 607 с. 2. Макаров Е.Г. Инженерные расчеты в Mathcad 14: учебный курс. СПб. : Питер, 2007. 592 с. 3. Трушин С.И. Метод конечных элементов. Теория и задачи: учебное пособие. М. : Изд-во АСВ, 2008. 256 с. 4. Хечумов Р.А. Применение метода конечных элементов к расчету конструкций / Р.А. Хечумов, Х. Кепплер, В.И. Прокопьев. М. : Изд-во АСВ, 1994. 353 с. 30 Учебное издание Игнатьев Александр Владимирович, Михайлова Наталия Анатольевна, Ерещенко Татьяна Владимировна МЕТОД КОНЕЧНЫХ ЭЛЕМЕНТОВ И ЕГО РЕАЛИЗАЦИЯ В СРЕДЕ MATHCAD Лабораторный практикум по дисциплине «Аналитические и численные методы решения уравнений математической физики» Зав. РИО О.Е. Горячева Зав. редакцией М.Л. Песчаная Редактор О.А. Шипунова Компьютерная правка и верстка Н.А. Дерина Подписано в печать 30.06.10. Формат 60х84/16. Бумага офсетная. Печать трафаретная. Гарнитура Таймс. Усл. печ. л. 1,9. Уч.-изд. л. 1,7. Тираж 100 экз. Заказ № 70 Государственное образовательное учреждение высшего профессионального образования «Волгоградский государственный архитектурно-строительный университет» Редакционно-издательский отдел Сектор оперативной полиграфии ЦИТ 400074, Волгоград, ул. Академическая, 1 31

Здесь приводится изложение упрощенного алгоритма решения плоской задачи механики деформируемого твердого тела методом конечных элементов в пакете Mathcad, опубликованного в моей статье в журнале Exponenta.Pro (№3, 2003 г.), а также на форуме Exponenta.ru. Дату поста ставлю соответствующую.

Рассматривается простейший и в то же время наиболее распространенный вариант с разбиением области на треугольные элементы. (По ходу дела ориентировался на алгоритм, изложенный в книге Фадеев А.Б. Метод конечных элементов в геомеханике. - М.: Недра, 1987 ).

1. Подготовка исходных данных.

Поскольку необходимо задать информацию о каждом из элементов и узлов расчетной области, удобнее всего использовать для подготовки данных табличный редактор Excel, тем более что в Mathcad предусмотрена возможность импорта данных из файлов формата.prn. Создается в Excel два файла с таблицами, содержащими сведения об элементах и узлах. Структура таблиц и размерности величин в них приведены на рис. 1 и 2. В таблице данных об узлах имеются два столбца специальных переменных Рх и Ру , которым присваивается признак фиксации перемещения вдоль оси 0х или 0у соответственно (принимает значение 1, если задано нулевое перемещение и 0 – если перемещение неизвестно).

Рис. 1. Структура таблицы исходных данных с информацией об элементах.

Рис. 2. Структура таблицы исходных данных с информацией об узлах и заданных узловых силах и перемещениях.

Для сохранения таблиц в нужном формате, выбираем Файл->Сохранить как… , указываем в соответствующих полях имя файла и тип – Форматированный текст (разделитель – пробел) . После нажатия кнопки Сохранить , нажать в появившемся окне Да . Таким образом, получаем файлы с именами, например, EL_1.prn и KN_1.prn.

2. Загрузка данных в Mathcad. Подготовка переменных.

Для удобства нумерации элементов массивов далее в книге Mathcad индекс первых элементов в массивах устанавливается равным единице:

Для получения данных из файлов в Mathcad используется функция READPRN(“filename.prn”) (возможно указание полного пути к файлу, в противном случае используется текущая папка, путь к которой можно узнать с помощью функции CWD).

Предположим, ранее созданные файлы находятся в папке DATA на диске D: Присваивается их содержимое матрицам DEL и DKN:

Присвоим соответствующим переменным значения из матриц:

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

Площадь n–го элемента удобно задать в виде функции пользователя (в векторе V перечисляются глобальные номера узлов элемента):

Общая площадь расчетной области:

В массовые силы пересчитывается вес элементов, поровну на каждый из их узлов:

Узловые силы, перемещения и их признаки размещаются в векторах последовательными парами значений: на четных позициях вертикальные компоненты, на нечетных горизонтальные:

3. Расчет матрицы жесткости системы.

Матрица жесткости системы получается путем объединения матриц жесткости элементов [K] , которые, в свою очередь рассчитываются по следующему выражению

Где Delta - площадь элемента; [B] - матрица производных функций формы (функция влияния узлов), [D] - матрица связи напряжений и деформаций:

Площадь элемента вычисляется ранее заданной функцией пользователя A(n) . Матрицу [D] также удобно задать в виде функции пользователя; для условий плоской деформации она будет иметь вид

Матрица [B] связывает между собой перемещения узлов элемента с его деформацией:

(выражения для функций формы Nj , Nk получаются путем круговой подстановки индексов в порядке i , j , k ),

i, j, k - номера узлов элемента, xi,j,k , yi,j,k - координаты узлов.

После несложных преобразований матрицу [B] можно представить в виде

Матрицу [B] представим в виде функции пользователя, задав предварительно вспомогательную матрицу P , определяющую порядок перестановки индексов в функциях формы:

Матрица жесткости системы вычисляется в следующем программном блоке:

(объединение матриц жесткости элементов в МЖС производится по следующему правилу: член МЖС Kci,j является суммой членов Ki,j из матриц жесткости всех элементов, примыкающих к узлу с і-й степенью свободы).

4. Решение системы уравнений

После этого і-й столбец и і-я строка МЖС, а также і-й неизвестный член в векторе сил могут быть удалены. Для удаления строк и столбцов из МЖС используем субматрицы, заданные функциями пользователя; М11 - удаляет первую строку и столбец, Mnn - последние, МI-IV - промежуточные.

Таким образом, нужно решить систему линейных алгебраических уравнений (СЛАУ) . В данном случае возможности системы Mathcad позволяют сильно упростить задачу. Для этого предусмотрена функция lsolve(M,V) для нахождения вектора решения СЛАУ, коэффициенты которой содержатся в массиве М, а свободные члены – в векторе V.

Программный модуль слева возвращает «на места» в общем векторе заданные узловые перемещения, ранее из него удалённые. Второй блок создаёт два вектора с осевыми компонентами узловых перемещений.

5. Нахождение осевых деформаций и напряжений в элементах

Зная полученные узловые перемещения, можно рассчитать для каждого элемента деформации и напряжения по соотношениям, упомянутым в п.3 (сигма и эпсилон):

В каждом элементе также подсчитываются главные напряжения и угол между осью 0у и вектором максимального главного напряжения . Чтобы избежать деления на ноль в строке вычисления угла использовано условное выражение, которое в случае равенства нулю выражения в знаменателе дроби присваивает углу значение .
.

6. Сохранение результатов.

Расчет по вышеизложенной процедуре занимает довольно непродолжительное время (например, на ПК с процессором Pentium-IV-1300 MHz; 128 MB RAM время расчета для области из 119 элементов 95 узлов составляет ~3 сек.), тем не менее, желательно сохранить результаты для последующего анализа.

Для этого сформируем матрицы, характеризующие НДС и поле перемещений, записав в них также координаты центров элементов и узлов:

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

Для записи данных в файл в Mathcad предусмотрена функция WRITEPRN(«filename.prn»); перед её использованием можно задать предварительно количество знаков после запятой переменной PRNPRECISION и ширину столбца в файле переменной PRNCOLWIDTH:

Рис. 3. Расчетная схема и её конечно-элементное представление.

В данном случае при разбиении на треугольные элементы получилась сеть из 95 узлов и 119 элементов. Нумерация – произвольная.

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

В силу симметрии граничные условия по перемещениям следующие: горизонтальные компоненты вдоль вертикальной (x=0 ) и вертикальные вдоль горизонтальной (y=0 ) сторон квадрата равны нулю. Неизвестны перемещения всех узловых точек внутри массива, на контуре выработки и на грани области.

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

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

Рис. 5. Примеры визуализации результатов расчета.

Литература:

1. Фадеев А.Б. Метод конечных элементов в геомеханике. – М.: Недра, 1987. – 221с.
2. Ержанов Ж.С., Каримбаев Т.Д. Метод конечных элементов в задачах механики горных пород. – Алма-Ата: Наука, 1975. – 239 с.
3. Зинкевич О. Метод конечных элементов в технике: Пер. с англ. – М.: Мир, 1975. - 542 с.
4. Норри Д., де Фриз Ж.. Введение в метод конечных элементов: Пер. с англ. – М.: Мир, 1981. – 304с.
5. Carlos A. Felippa. Introduction to Finite Element Methods. – Department of Aerospace Engineering Sciences and Center for Aerospace Structures University of Colorado, Boulder. – 2001.
6. Kyran D. Mish, Leonard R. Herrmann, LaDawn Haws. Finite Element Procedures in Applied Mechanics (попадалось где-то в инете).
7. Зенкевич О., Морган К. Конечные элементы и аппроксимация. – М.: Мир, 1986. – 318 с.
8. Зенкевич О., Чанг И. Метод конечных элементов теории сооружений и в механике сплошных сред. – М.: Недра, 1974. – 240 с.

Ссылки:

• http://www.fea.ru/ ...Сайт FEA.RU посвящен актуальным проблемам конечно-элементной механики и компьютерного инжиниринга (CAE), МКЭ и расчётам на прочность;
• http://www.cae.ru/ форум CAD и CAE-систем в том числе теоретические и прикладные аспекты КЭ моделирования и решения задач механики деформируемого твердого тела. Механика конструкций, машин, сооружений и установок;
• - Мощный каталог ресурсов, касающихся МКЭ;
• http://www.isib.cnr.it/~secchi/EdMultifield/ - Сайт программы для конечноэлементных расчетов с неплохим описанием метода.