CC BY
132
Если бы мы в соответствии с рекомендациями, данными в [2], стали проводить эту работу для достаточно большого числа гармоник ряда (1), отмечая влияние саморегулирования для каждой из них, то нет уверенности в том, что, просуммировав полученные решения, мы пришли бы к установленному результату, полученному прямым интегрированием уравнения (2). Так что при использовании в расчетах статической характеристики электродвигателя саморегулирование проявляется не для всякого вида нагрузки. Не исключено, впрочем, что оно окажется заметным при применении динамической характеристики электродвигателя [4].
Отмечена также тенденция к увеличению коэффициента неравномерности движения при уменьшении времени действия момента полезного сопротивления (или угла (рм) как при учете саморегулирования, так и без него даже при неизменном импульсе этого момента 5м = -хТМ (или работе
Ам = -Мрм&-200ж Дж). Следуя [1], можно таким же образом улучшить ряд Фурье для производной функции, если она окажется кусочно-
непрерывной: / (г) = (г) = /2 (г) + ф (г), где
ф(г) - вспомогательная кусочно-постоянная функция. Тогда после интегрирования получится представление для исходной функции
/ (г) = ФМ)+¥г(г) + ),
где ^(г) = |/2(г )&,, Ф(г) = \ф(г )& - первооб-
разные. Оно включает ряд Фурье для Г2(г) с коэффициентами, имеющими порядок убывания не ниже —N . Постоянная интегрирования является нулевым его членом и находится обычным способом:
С = ^ = 2
1 т 1 т
-1Ф1 (г )Л = - (г) - ф0 (г) - ^ (г )]л.
т о т о
Правда, ф (г) будет уже кусочно-линейной функцией. Так что при последовательном применении этого метода вместе с улучшением сходимости тригонометрического ряда усложняются вспомогательные функции, что может вызвать затруднения с интегрированием дифференциальных уравнений.
Таким образом, если не обеспечена малость влияния высших гармоник ряда Фурье, в который разложено входное (внешнее) воздействие на механическую систему: оно предполагается близким к синусоидальному или передаточная функция системы «гасит» высшие гармоники [5], - то надо внимательно отнестись к характеру этого воздействия. Если оно представляет собой разрывную периодическую функцию, то необходимо улучшать сходимость ее ряда Фурье описанным выше способом, иначе результаты, полученные при использовании одной или нескольких первых гармоник ряда, могут оказаться ошибочными.
Библиографический список
1. Смирнов В.И. Курс высшей математики: учеб. для вузов. М.: Наука, 1974. Т. 2. 656 с.
2. Фролов К.В. Теория механизмов и машин: учеб. для втузов. М.: Высшая школа, 1987. 367 с.
3. Клепацкий А.Н. Оценка влияния саморегулирования на неравномерность движения машин с кратковременно-повторным действием сил полезного сопротивления // Акту-
альные вопросы прикладной науки и образования: сборник. Иркутск: Изд-во ИрГТУ, 2010. С. 33-40.
4. Коловский М.З. Динамика машин. Л.: Машиностроение, 1989. 263 с.
5. Попов Е.П., Пальтов И.П. Приближенные методы исследования нелинейных автоматических систем. М.: Физ-матгиз, 1960. 791 с.
УДК 621.81
РАЗРАБОТКА ПРОГРАММЫ РЕАЛИЗАЦИИ МЕТОДА КОНЕЧНЫХ ЭЛЕМЕНТОВ В ЗАДАЧЕ СТАЦИОНАРНОЙ ТЕПЛОПРОВОДНОСТИ
А
А.А.Кудрявцев1
Национальный исследовательский Иркутский государственный технический университет, 664074, г. Иркутск, ул. Лермонтова, 83.
Рассмотрена задача трехмерной стационарной теплопроводности с решением на основе метода конечных элементов (МКЭ). Представлена блок-схема алгоритма решения задачи стационарной теплопроводности на основе МКЭ с программированием на языке Фортран. Приведено сравнение результатов с теоретическим решением. Предложено разработать комплексную математическую модель для решения контактной задачи теплонапряжен-ности сборных конструкций. Ил. 3. Табл. 1. Библиогр. 5 назв.
Ключевые слова: теплопроводность; метод конечных элементов; программирование.
1-
Кудрявцев Александр Александрович, старший преподаватель кафедры самолетостроения и эксплуатации авиационной техники, тел.: (3952) 405133, 89021786789.
Kudryavtsev Alexander, Senior Lecturer of the Department of Aircraft Construction and Operation of Aircrafts, tel.: (3952) 405133, 89021786789.
DEVELOPMENT OF THE FINITE ELEMENT METHOD REALIZATION PROGRAM IN THE PROBLEM OF STATIONARY HEAT CONDUCTIVITY A.A. Kudryavtsev
National Research Irkutsk State Technical University, 83, Lermontov St., Irkutsk, 664074.
The problem of three-dimensional stationary heat conductivity with the solution based on the finite element method (FEM) is considered. A block diagram of the algorithm to solve the problem of stationary heat conductivity on the basis of FEM with the Fortran programming is presented. The results are compared with the theoretical solution. It is proposed to develop a comprehensive mathematical model for solving the contact problem of thermal stress of prefabricated structures.
3 figures. 1 table. 5 sources.
Key words: heat conductivity; finite element method; programming.
Актуальность анализа теплонапряженности сборных конструкций турбомашин обусловлена достижением в них высоких уровней температур, напряжений и деформаций в процессе работы. Полное решение представленной проблемы требует проведения комплекса взаимозависимых исследований, включая решение задач термодинамики и теории упругости. Наряду с этим, решение задач стационарной теплопроводности является неотъемлемой частью этих исследований, а повышение точности их решения определяет актуальность представленной ниже работы.
Общее уравнение теплопроводности имеет вид [4]:
/-\2тг
т^ O T O T „ O T „ .
Kxx + Kyy + Kzz + Q = 0 , (1)
Ox Oy Oz
где T- температура, °С (°К); Kxx, Kyy, Kzz - коэффициенты теплопроводности в направлениях x, y и z, Вт/м ■ K; x, y, z - геометрические координаты декартовой системы координат; Q - внутренний источник тепла (при наличии), Вт/м3.
Если на границе тела происходит конвективный теплообмен или задан поток тепла, то с уравнением (1) связывают граничное условие
Kxx ~rlx + Kyy ^ly + Kzz + h(T - T ) + q = 0, (2)
Ox Oy Oz
где h - коэффициент теплообмена, Вт/м2-К; Г» - температура окружающей среды, К; lxx, lyy, lzz - направляющие косинусы вектора нормали к поверхности; q - удельный тепловой поток (плотность теплового потока), Вт/м2, который считается положительным, если тепло теряется телом.
Для тел сложной формы уравнение теплопроводности практически невозможно решить аналитически, поэтому применяют численные методы, среди которых одним из наиболее мощных является метод конечных элементов (МКЭ). Этот метод в данной работе взят за основу.
МКЭ основан на аппроксимации сплошной среды с бесконечным числом степеней свободы совокупностью подобластей (конечных элементов), имеющих конечное число степеней свободы. В каждой из этих подобластей распределение неизвестной функции аппроксимируется базисными функциями, построенными с использованием простейших алгебраических полиномов и выраженными через значения аппроксимируемой функции в узлах конечного элемента. Алгебраические полиномы, используемые при построении функций форм, выбираются с учетом условия непрерывности решения в межэлементной зоне. Это необходимо для достижения сходимости приводимого численного решения к точному решению.
Дальнейшее решение задачи МКЭ сводится к вариационной минимизации функционала рассматриваемой задачи. Это позволяет получить приближенное решение, не прибегая к интегрированию системы дифференциальных уравнений.
С вариационной точки зрения решение уравнения (1) эквивалентно отысканию минимума функционала [3]:
' = •2
V 2
i 1
Kx If J +
2 , Лг
+
Oy
> qT +1 h(T - Tj
OT + Kzz I OH- 2QT
dV +
(3)
dS.
5 - ^
Минимизация функционала (3) при реализации МКЭ осуществляется на множестве узловых значений {7} конечно-элементной модели. Процесс минимизации начинается с преобразования функционала путем введения двух матриц. Члены первой матрицы представляют собой величины первых производных неизвестного поля температур относительно координатно-геометрических характеристик в декартовой системе координат. Общее выражение для этой матрицы имеет вид
ыт=
(4)
дТ дТ дТ дх ду дz
Члены второй матрицы содержат коэффициенты теплопроводности материала в используемой трехосной
системе координат. Она называется матрицей свойств теплопроводности и имеет следующий вид:
IX 0 0 -
М = 0 Куу 0 . _ 0 0 К„ _
лено как
■Н\
Теперь соотношение (3) может быть представлено как
X = 11Ы№}-ЩТ^У + |qTdS + 11[т2 - 2ТТШ + Т2.
$1 $2
(5)
(6)
у ~
Здесь следует отметить, что поток тепла ц и конвективная потеря тепла Л(Т-Т„) не имеют места на одном и том же участке поверхности границы. Если существуют потери тепла за счет конвекции, то отсутствует отвод или приток тепла за счет теплового потока и наоборот [4].
Так как при использовании МКЭ функция температуры Т не является непрерывной по всей области определения, необходимо ввести в рассмотрение функции Т;е), определенные и непрерывные на отдельных конечных элементах (е). Это позволяет записать выражение вида
Х =
]Г | ^ } [Б{е) ]{$(е) }у - | а{е)Т{е^У + | ^е)Т(е^ +
е=1 ,Ле) 2 -<'*>
+
. и(е) г 1
Г-[т(е)Т(е) - 2Т(е)Тх + Тш2
,(е) 2
(7)
(е) $2
где Е - общее число конечных элементов. Выражение (7) может быть записано в виде
х = х(1)+х(2) +...+х( Е) = £х(е),
(8)
где Xе - вклад отдельного элемента в общий температурно-энергетический баланс рассматриваемого объекта, характеризуемый функционалом х
Минимизация функционала хтребует выполнения соотношения
дх д Л (е) £ дх"'
= 0.
(9)
д{т} д{Т£ д{т}
Чтобы определить частные производные в выражении (9), необходимо интегралы, входящие в уравнение (7), выразить через узловые значения температуры {Т}.
Интерполяционный полином, аппроксимирующий температуру в области определения конечного элемента, представляется в виде
X
Т(е) = [ы(е) ]{Т }=[ы(е) ы2е) ы2е) • •• ы(ре) ]
где Ыр - функция формы элемента в узле р; р - число узлов элемента. Тогда выражение для {д(е)} можно записать в виде
Ь<" }=
дТ(е) дЫ2е) ды(е дЫ(е)
дх дТ(е) дх дЫ(е) дх дЫ(2е) дх дЫ(е)
ду дТ(е) ду дЫ(е) ду дЫ(е) ду дЫ(е)
[ дz } дz дz дz
Т
Т
(10)
(11)
Е
е
2
>
3
2
или в сокращенном виде
{?(е) }=[в(е) ]{Т }, (12)
где [Б(е>] - матрица, содержащая информацию о первых производных функций формы в глобальной системе координат.
Поскольку функции формы задаются в локальных координатах, необходимо выразить глобальные производные, входящие в (11), через локальные производные [1]. Это осуществляется через матрицу Якоби
дх ду &
дх ду дг
дц дц дц
дх ду дг
дС дС ~дС_
где £ п, С - локальные координаты конечного элемента.
Таким образом, глобальные производные функций формы можно представить в виде
[/] =
(13)
'дЫ' 'дЫ'
дх
дЫ дЫ
ду дц
дЫ дЫ
1& \ [дС \
(14)
Выражение (10) можно переписать в следующем виде:
Т(е) = [ы(е) ]{Т }={Т }Т [ы(е) ]. Формулы (10), (11) и (12) позволяют выразить интегралы по элементам в (7) в виде
,«- | 1 {Т}Т[в(е)][&е)][в(е)]{Т\1У - |о[м(е)]{Т}у
(15)
^е) =
+
+
|д[м(е) ]{т}5 + | Н {т}Т[Ы(е)] ]{Т}5
(16)
5
| нф(-е) ]№ + | Н ^.
Проводя дифференцирование функционала (16) относительно величины температуры на отдельном конечном элементе, рассмотрим выражения для каждого из его составляющих:
| 1{т}Т[в(е)][о(е)\в(е)\T\iV = |[в(е)][о(е)\в(е)\T\iV;
д{Т}у(е) 2 V<е)
-д-1 ^ = | еИ ]dv;
д{Т }Vе V(е)
-д. \ ^ = \ ^ ]dS ;
-д- \ н{Т}ТИ^ =\н[Ые)ГИ^ ;
д{Т } 5« 2 5«
-¿т |нТ.И= |НТN1[dS;
д{Т } 5(е 5(е)
[ = 0. дТ} 1 2 -
(17)
<
(е)
5
Подставляя выражение каждого из составляющих (17) в выражение для производной дХе/д{Т на отдельном конечном элементе, имеем соотношение
д{Т }
| [д(0 ] [D(e) Iв(e) + | к^(e) ] [N(e) №
у V(e > ¿(е >
{T }-
. (18)
- | б^(е)^ + | q[N(e) ]dS - |кТт [N(e) ]dS.
V
(е) + \q\NieЦdS - \кТх ^^(е)
(е) о(е) о(е)
Эта совокупность интегралов может быть записана в компактной форме
д^(е)
дТ}
= [к(е) ]{т}+{/ (е)}, (17)
где
[к(е) ] = |[я(е) ][^е) ][я(е) ¡V + | к[^е) ] [N(e) ^ (20)
V
(е) х(е)
{¡(е) }=-{ XdV + | q[N(e) ]dS - | кТ^^ ] dS. (21)
5 (е) S (е)
Интегралы в выражениях (20) и (21) определяют соответственно матрицу теплопроводности элемента [Л(е)] и
вектор нагрузки (теплообмена) элемента {/е)}. Для вычисления этих интегралов необходимо использовать приемы численного интегрирования [1, 4]. При этом используется квадратурная формула Гаусса-Лежандра. Для интегрирования выражений (20) и (21) достаточно квадратуры второго порядка [4].
Глобальная система алгебраических уравнений МКЭ относительно температурного состояния рассматриваемого объекта моделирования получается после подстановки выражения (19) в выражение (9):
-£(к«]т}+{/")8=о,
д{г}
или
[Г]{Т}={^}, (22)
где [К ] = ^[к(е) ] - глобальная матрица блочно-диагонального типа, отражающая свойства теплопроводности
е=1
отдельных составляющих сборную конструкцию твердых тел; {Т - глобальный вектор-столбец неизвестных значений температуры рассматриваемого объекта, получаемый после решения задачи стационарной теплопроводности и теплообмена; {^ }= {/(е)} - глобальный вектор-столбец теплообмена с внешней средой.
е=1
Построение глобальной матрицы теплопроводности, то есть формирование ее суммированием из отдельных матриц теплопроводности элементов, осуществляется аналогично методу прямой жесткости при решении силовой задачи. Подобная процедура применяется и при построении вектора теплообмена.
Решением глобальной системы алгебраических уравнений (22) имеем поле температур, представленное в узловых точках значениями коэффициентов глобального вектора {Т}.
Представленная математическая модель на основе МКЭ была запрограммирована на алгоритмическом языке ФОРТРАН, возможности которого вполне достаточны для программирования при решении задач МКЭ. Алгоритм программы представлен на рис. 1.
Исходные данные для ФОРТРАН-программы считываются в автоматическом режиме с электронной модели, построенной в программе подготовки конечноэлементных моделей РБМДР. Файл результатов, создаваемый ФОРТРАН-программой после расчетов, можно также открывать в программе РБМДР для более удобного и наглядного просмотра. Данная реализация решения задачи стационарной теплопроводности с программированием на языке ФОРТРАН необходима для дальнейшего включения в этот алгоритм математического аппарата по решению контактной задачи теплонапряженности сборных конструкций. Это позволит создать комплексную математическую модель, включающую математические модели теплопроводности и напряженно -деформированного состояния деталей сборных конструкций. Такая комплексная модель является необходимым элементом системы инженерного анализа высоконагруженных тепловых машин.
В качестве тестового примера анализа трехмерной стационарной теплопроводности рассматривается куб (рис. 2), грани которого имеют постоянную температуру, причем одна грань имеет температуру 100 °С, а пять других - 0 °С. Рассматриваемая задача уже имеет численное (методом конечных разностей) и аналитическое решения [5], причем аналитическое решение для данной задачи имеет вид бесконечного сходящегося ряда [2]:
и
( е)
1 .Ввод исходных данных с файла *.(М
1 г
2.Вычисление ненулевых членов матрицы [О]
1 г
З.Задание координат точек интегрирования
1 г
4.Цикл по конечным элементам
б.Цикл по точкам интегрирования
1
6.Вычисление производных функций формы по координатам точек интегрирования
1
7.Вычисление матрицы Якоби
1 г
8.Вычисление определителя матрицы Якоби
Э.Вычисление обратной матрицы Якоби
10.Вычисление производных функций формы по глобальным координатам
11 .Вычисление элементарного объема
12.Вычисление ненулевых членов матрицы [6]
13.Транспонирование матрицы [В]
14.Вычисление матрицы теплопроводности конечного элемента
15.Суммирование элементов матриц [Г>] с учетом глобальной нумерации узлов -построение матрицы [К]
16.Формирование вектора тепловых нагрузок {/=}
I
17.Учет граничных условий в матрицах [К] и {У1}
I
18.Модификация матриц [К] и {/=} по Гауссу
19.Решение системы уравнений - нахождение вектора {7}
20.Формирование файла результатов *.ГО6
Рис. 1. Блок-схема алгоритма решения задачи стационарной теплопроводности на основе МКЭ
16 " " 100вЬ[/(а - х)]б1п[(2р + 1)лу / а]вт[(2р + 1)ж / а]
ж
р=0 9=0
(2 р + 1)(2д + 1>Ь(/а)
72
где / = ратура 7.
(2 р + 1)2ж2 (29 + 1)2ж2 ---^--ь--^-; а - ребро куба; х, у, г - координаты точки, для которой находится темпе-
а
а
Рис. 2. Пространственная сетка куба (представлена % куба ввиду симметричности); наклонным шрифтом показаны граничные условия, прямым шрифтом - номера узлов, в которых определялась температура
Результаты решения представленной задачи стационарной теплопроводности, полученные по ФОРТРАН -программе, показаны на рис. 3 в виде изолиний для % куба ввиду симметрии температурного поля.
Рис. 3. Распределение температур в % куба, °С
Значения температур для узлов сетки, пронумерованных, как показано на рис.2, сведены в таблицу. Решения представлены по методу конечных разностей [5], аналитически [2], а также по МКЭ, причем для МКЭ приводится несколько вариантов в зависимости от густоты сетки. _Значения температур в узлах сетки, °С_
Номер узла Метод конечных разностей Аналитический метод МКЭ
4 элемента на ребре Отклонение относительно аналитического метода, % 8 элементов на ребре Отклонение относительно аналитического решения, % 16 элементов на ребре Отклонение относительно аналитического решения, %
1 30,4 30,7 30,4 0,98 30,8 0,33 30,7 0,0
2 9,8 9,1 9,2 1,1 8,8 3,3 9,0 1,1
3 3,0 2,6 2,9 11,5 2,5 3,8 2,6 0,0
4 36,2 37,3 37,2 0,27 37,6 0,80 37,4 0,27
5 12,8 12,3 11,5 6,5 12,1 1,63 12,2 0,81
6 4,0 3,7 3,6 2,7 3,6 2,7 3,6 2,7
7 36,2 37,3 37,2 0,27 37,6 0,80 37,4 0,27
8 12,8 12,3 11,5 6,5 12,1 1,63 12,2 0,81
9 4,0 3,7 3,6 2,7 3,6 2,7 3,6 2,7
10 43,6 45,8 49,1 7,2 46,6 1,75 46,0 0,44
11 16,7 16,7 16,7 0,0 16,7 0,0 16,7 0,0
12 5,5 5,1 4,8 5,9 5,0 1,96 5,1 0,0
Из приведенных решений можно констатировать совпадение результатов, полученных реализацией алгоритма, представленного выше, с отклонением, не превышающим 11,5% по отношению к аналитическим данным, даже при достаточно крупном разбиении (4 элемента на ребре). С увеличением густоты сетки точность в сравнении с аналитическими данными повышается, однако использование сетки с количеством элементов на ребре куба более 16 уже нецелесообразно, так как выигрыш точности не оправдывается потерями времени на расчеты.
Библиографический список
1. Зенкевич О.С. Метод конечных элементов в технике. М.: Мир, 1975. 542 с.
2. Карслоу Г., Егер Д. Теплопроводность твердых тел. М., 1964.
3. Пыхалов А.А., Милов А.Е. Контактная задача статического и динамического анализа сборных роторов турбомашин: монография. Иркутск: Изд-во ИрГТУ, 2007. 192 с.
4. Сегерлинд Л. Применение метода конечных элементов. М.: Мир, 1979. 392 с.
5. Юдаев Б.Н. Теплопередача: учебник для вузов. 2-е изд., перераб. и доп. М.: Высш. школа, 1981. 319 с. УДК 620.17
РЕГРЕССИОННЫЙ АНАЛИЗ СИЛЫ УДАРНОГО ВЗАИМОДЕЙСТВИЯ УПРУГО-ВЯЗКО-ПЛАСТИЧНОЙ МЕХАНОРЕОЛОГИЧЕСКОЙ МОДЕЛИ
В.Л.Лапшин1, А.В.Глухов2
Национальный исследовательский Иркутский государственный технический университет, 664074, г. Иркутск, ул. Лермонтова, 83.
Рассматривается упруго-вязко-пластичная механореологическая модель, предназначенная для исследования процессов ударного взаимодействия тел. Модель позволяет изучить влияние упругих, вязких (диссипативных) и пластических свойств материала на динамику ударного взаимодействия тел. Приводятся результаты компьютерного исследования, отражающие влияние упругих и пластических параметров модели на силу ударного взаимодействия, полученные путем проведения факторных экспериментов с использованием уравнений регрессии. Ил. 8. Табл. 3. Библиогр. 9 назв.
Ключевые слова: упруго-вязко-пластичная модель; механореологическая модель; ударное взаимодействие тел; компьютерное моделирование.
1-
Лапшин Владимир Леонардович, доктор технических наук, профессор, зав. кафедрой сопротивления материалов и строительной механики, тел.: (3952) 405425, е-mail: lapshin@istu.edu
Lapshin Vladimir, Doctor of technical sciences, Professor, Head of the Department of Strength of Materials and Structural Mechanics, tel.: (3952) 405425, e-mail: lapshin@istu.edu
2Глухов Александр Владимирович, студент, тел.: (3952) 405425, е-mail: lapshin@istu.edu Glukhov Alexander, Student, tel.: (3952) 405425, e-mail: lapshin@istu.edu
