Математическое моделирование термоупругого контактного взаимодействия осесимметричных тел Текст научной статьи по специальности «Физика»

УДК 536.2+539.3

Математическое моделирование термоупругого контактного взаимодействия осесимметричных тел

© С.М. Богатырь1, М.П. Галанин2'3, В.И. Кузнецов1, В.В. Новиков1, А.С. Родин2, М.Е. Яковлев3, А.В. Крупкин1, В.В. Лукин2'3, И.В. Станкевич3

1 ОАО «ВНИИНМ им. акад. А.А. Бочвара», Москва, 123098, Россия 2 Институт прикладной математики им. М.В. Келдыша РАН, Москва, 125047, Россия 3 МГТУ им. Н.Э. Баумана, Москва, 105005, Россия

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

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

Введение. Для надежной оценки ресурса элементов конструкций объектов энергомашиностроения, работающих в условиях высокоинтенсивного термомеханического нагружения, вызванного в том числе и контактным взаимодействием, важным является оценка напряженно-деформированного состояния. Таким образом, возникает необходимость решения контактных задач механики деформируемого твердого тела (МДТТ) с учетом особенностей их конструкции и условий термомеханического нагружения. Аналитические решения контактных задач получены для весьма ограниченного числа видов контактного взаимодействия и форм контактирующих поверхностей. В большинстве практически важных ситуаций, связанных с принятием конструктивных решений, например, для контактирующих тел, имеющих сложную геометрическую форму, и при сравнительно невысоких требованиях к гладкости функций, входящих в формулировку краевых и начально-краевых задач, наиболее перспективны численные методы. В настоящее время для численного решения контактных задач широко применяется конечно-элементная технология, в рамках которой реализуется выбранный алгоритм. Основными алгоритмами решения контактных задач являются метод множителей Лагранжа, метод штрафов и их комбинации, а также релаксационные схемы [1].

При решении некоторых типов задач в рассмотрение вводят «псевдосреду» [2]. Весьма перспективным для решения контактных задач является применение альтернирующего метода Шварца, основанного на принципе поочередности [3-7]. Преимущества этого метода состоят в том, что не требуется согласовывать построение узлов конечно-элементных моделей на поверхностях контакта и переформировывать матрицы систем линейных алгебраических уравнений в процессе итерационного уточнения границ зон контакта. Ниже рассмотрены некоторые особенности алгоритмического оформления альтернирующего метода Шварца применительно к решению осесимметричных контактных задач теории термоупругости, имеющих выраженную практическую направленность. Необходимо отметить, что алгоритмы, используемые для решения контактных задач теории термоупругости, служат основой для решения контактных задач механики деформируемого твердого тела с учетом деформаций пластичности и ползучести.

Математическая постановка осесимметричной краевой задачи термоупругости с учетом контактного взаимодействия. Рассмотрим два осесимметричных линейно-упругих контактирующих тела А и В, находящиеся в условиях о се симметричного термомеханического нагружения. Если зафиксировать некоторую плоскость я, в которой лежит ось вращения I, то исследуемые тела геометрически можно охарактеризовать двумерными областями Са, а Е [А, В}, лежащими в плоскости я и ограниченными кусочно-гладкими границами дСа. Введем цилиндрическую систему координат Огфг, причем ось Ох совпадает с осью вращения I.

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

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

Рассмотрим следующую начально-краевую задачу для нелинейного уравнения теплопроводности:

дТ

сЕ(М,Т)р(М,Т) — = {кгз (М,Т)Т,3) п +ду (М,г,Т), (1)

(М,г) Е йа х (0,гм];

Т(M, 0) = То (M), M G Ga; (2)

т(M,t)\sa = Tw(M,t), M G sa С dGa, t> 0; (3) -пгкгз(M,T)T,3\Saq = qw(M,t,T), M G s' С dGa, t> 0; (4)

-пгкгз(M,T)T,j \Sq = aw(M,t,T)(T(M,t) - Tf (M,t)),

Tf

M G sqf С dGa, t > 0, (5)

где t — время; T0(M) — начальная температура; T(M,t) — температура в момент времени t; kij (M,T) — компоненты тензора теплопроводности; qv(M,t) — мощность внутренних источников (стоков); ce(M,T) — удельная теплоемкость среды, занимающей область Ga; р(М, Т) — плотность среды, занимающей область Ga; Tw(M, t) — температура поверхности Sa; qw (M,t,T) — плотность теплового потока на поверхности S^; aw(M,t,T) — коэффициент теплоотдачи на поверхности Sa ; Tf (M, t) — температура среды у поверхности STf ; щ — компоненты единичного вектора внешней нормали n = niei к границе dGa. Здесь и ниже a G [А, В}, i,j = r,z, а запятой с индексом обозначена операция дифференцирования по пространственным координатам Г, Z.

Предположим, что между контактными поверхностями двух тел a, b G [А, В}, a = b, S'a = Sf = Sk имеется неидеальный тепловой контакт, тогда

щЬ' (M, Ta) Ta,j |s = afc (M) (Tb (M, t) - Ta (M, t)) ; (6) rnka (M, Ta) Ta,j | sfc = Щkl (M, Tb) Tb,3 | sfc . (7)

Здесь ni — компоненты вектора внешней нормали к контактной поверхности Sa в точке M ; a^ (M ) — контактный коэффициент теплоотдачи; Ta (M,t), Tb (M,t) — температуры сходственных точек тел a и b на контактной поверхности, M G S к = Sa = Sf.

Отметим, что контактный коэффициент теплоотдачи a к (М ) может зависеть от контактного давления an, шероховатости поверхностей контакта и от других параметров [8].

Краевая квазистатическая задача МДТТ в данном случае формулируется так [9, 10]: уравнения равновесия

Sji,j [u (M) ,T (M,t)} + Qi (M) = 0, M G Ga, (8) граничные условия (кинематические и силовые соответственно)

u (M)\Si = u0 (M), M G s' С dGa, (9)

a^ (u, t) Uj\Sa = Pi (M), M g s; С dGa, (10)

соотношения Коши

ди дь и ди

I ди ди

е(М ) = {е} = К , в, , £ф ,уГ2 }т = ^ — ,— ,- + , М е С0

[ ОТ ОХ Г ОХ ОТ )

(11)

и определяющие уравнения — закон Гука в виде

а(М) = {о>, ох, оф, огх}т = Н(е - е0) =

= [Н] ({ег, в,, вф, уГ2}т - {е0, в0, в£, у0*}Т) , М е Са, (12)

где Н = [Н] — матрица Гука; а = {ог, ох, о3, огх }т — вектор напряжений; е = {в} = {вг, вх, в3, угг}т — вектор деформации; е0 = {в0} = = {в0, е°, в3,у^}т — вектор начальной деформации (в данном слу-

3 ">

чае — температурной); Qi — компоненты массовой (объемной) силы Q (М) = = {дг}т; и (М) = {и(М),ь(М)}т — вектор перемещений точек тела, и (М) — перемещение в направлении оси Ог, V (М) — перемещение в направлении оси Ог; ы0 = {и0, ^0}т — вектор заданных перемещений точек поверхности Pí (х) — компоненты заданной распределенной нагрузки р = {р} = {рг ,рг}т на поверхности Б*.

Запись Оу {ы,Т} означает, что выполнены все необходимые преобразования, позволяющие тензор напряжений а выразить через вектор перемещения ы (М) и температуру Т (М,Ь). В дальнейшем аргументы у функций опущены, если это не приводит к противоречию. Кроме массовых Q (М) и поверхностных р (М) сил в отдельных пд точках внешних поверхностей дСа, а е {А, В}, с фиксированными координатами хк (к = 1, п%) могут быть заданы дискретные силы Я*.

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

и* (М) - < (М) = Ьп (М), М е Бк; (13)

и по напряжениям (силовое условие)

о^ (М) = -о£ (М) ^ 0, М е вк; (14)

о^ (М) = о? (М) = 0, М е Бк. (15)

Здесь м^, и% — проекции векторов перемещений граничных точек на направление внешней нормали п к границе тела А; 8п — начальное расстояние (зазор) по нормали между граничными точками тел А и В;

Ли ^ А И

О-, о^ — проекции векторов напряжений аЛ и а° на направление внешней нормали п к границе тела А; о^, о^ — проекции векторов напряжений аА и ав в направлении общей касательной х к контактной поверхности . Условие (15) указывает на отсутствие трения в зоне контактного взаимодействия.

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

Основные матричные соотношения метода конечных элементов. Задача теплопроводности. Для решения нелинейных нестационарных задач теплопроводности контактирующих тел (1)-(7) применим метод последовательных приближений [11], в соответствии с которым их решение заменяют решением серии линейных задач с коррекцией всех данных, зависящих от температуры. Для линеаризации зафиксируем в каждом теле а е {А, В} некоторую температуру Т*(М) (М е Са) (в общем случае в каждом теле свою) как параметр в коэффициентах уравнений (1) св (М,Т*), р (М,Т*), к^ (М,Т*), члене цу (М,Ь,Т*), учитывающем мощность внутренних источников (стоков) теплоты, в граничных условиях (3) и (4): кц (М,Т*), (М,Т*), а,ш (М,Т*) и в условиях на контактной поверхности (6) и (7): Щ (М,Т*). Отметим, что при формировании контактных граничных условий (6) для тела а необходимо использовать зафиксированные значения температур сходственных точек Тр* (М) тела р, т. е. принять Тр (М, г) = Тр* (М).

После решения линеаризованных задач для тел а е {А, В} функции Т* (М) корректируются: Т* (М) == Т (М) и, если сходимость не достигнута, вновь выполняется решение задач (1)-(7).

Численное решение каждой линейной задачи теплопроводности предполагает построение соответствующего дискретного аналога. При решении нестационарных задач теплопроводности конечно-элементную дискретизацию по пространству целесообразно выполнять на основе метода взвешенных невязок в форме Галеркина [11, 12]. Тогда для каждого контактирующего тела а е {А,В} линеаризованная нестационарная задача сводится к решению задачи Коши для линейного матричного дифференциального уравнения первого порядка

[С] {Т} + [К] {Т} = {Д} (16)

с начальным условием

{Т }и = №},

где {Т0} — проекция функции Т0 (М) (начального условия (2)) на узлы сетки конечно-элементной модели.

При построении уравнения (16) использованы интерполяционные

соотношения [12]:

'Т =[мУ]1а¥] {Т} , Т Ч^Па?] {Т} ,

{Т,г , Т,г }т = , ,г ]т[а()]{Т} = [ВЩа$]{Т}.

Здесь [М-у[^9)] — матрицы-строки функций формы соответственно для объемного и поверхностного интегрирования, используемые при решении температурных задач; [а—"], [аз"] — матрицы геометрических связей соответственно объемных и поверхностных (и контактных поверхностных) конечных элементов, используемые при решении температурных задач; [В(е)] — матрица градиентов конечного элемента; (е), (д) — идентификационные метки соответственно объемных и поверхностных (в том числе и контактных) конечных элементов; {Т} — глобальный вектор узловых температур тела а € {А, В}.

Глобальные матрицы теплоемкости [С] и теплопроводности [К], а также глобальный вектор тепловой нагрузки {Д} строятся для каждого контактирующего тела а € {А, В} по следующим формулам:

ке1

[С] = У^^у^Г^2*1 [<$],

е=1

ке1 / С \

К] = ^1ау]т( 2* [В (е)]т[Д(е) ][В(е) ]г ¿г ¿г)[а(у;)] +

+ £[4е)]т(2* [ а^П^М)[4е)] +

е=1 ^ 7-, '

г(е)

+ £[49)]т(2*1 а? [^П^М)[49)],

9=1 (е)

4е)

ке1

{Д} = ^[ау]т2* I <1у[Му]тЫЫг

А(е)

[

е=1

£ [4е)]т2*(/№[М^]тг ¿1- I а ^[М^Гг^У

ау) ь:

(д) >

£[49)]т2* Г а^Т^ГМ,

9=1 Ь^) ьк

/

где ке1, кдг, кк — число соответственно объемных, поверхностных и контактных поверхностных конечных элементов рассматриваемого тела а е [А, В}; [Д(е)] = 1(0е) ^ — матрица коэффициентов

теплопроводности (предполагается, что оси Ог и Ог — главные оси анизотропии); ¿(е)

— область интегрирования объемного конечного элемента, лежащая в плоскости я; — длина (в общем случае криволинейного) ребра, характеризующая поверхностный конечный элемент (д), принадлежащий поверхности Бд; Ь^ — длина ребра, характеризующая поверхностный конечный элемент, принадлежащий поверхности ; Ь^ — длина ребра, характеризующая контактный поверхностный конечный элемент, принадлежащий поверхности Бк; Т^д) — температура сходственной точки, расположенной на контактной поверхности соседнего тела р е [А, В}, то есть тела, с которым контактирует рассматриваемое тело а, а = р.

Контактная задача теории упругости. Решение этой задачи выполняется итерационно с помощью процедур альтернирующего метода Шварца. На каждой итерации решается специальным образом линеаризованная задача теории упругости. Используя вариационную формулировку линейных задач теории упругости и построение дискретизации по пространству в рамках конечно-элементной технологии, получаем следующее матричное уравнение относительно глобального вектора узловых перемещений [и}:

[К] [и} = [Я},

где глобальные матрица жесткости [К] и вектор нагрузки [Я} [9]:

ке1 л

[К] = ^[ау]т2я [В (е)]т[Я(е)][Б (е)]г ¿г

е=1 М

[К} [а^]т2ж^ у [Б(е)]т[Я(е)] [в0е) }Ыг(1г++

е=1 А(е)

каг кд

+ ^ 1^д)]т{р(д)} гШ + ^ 2яг(1)[Ъ(1)]т{П(1)} .

д=1 А) =

Здесь [Ж^е)], [^д)] — матрицы-строки функций формы соответственно для объемного и поверхностного интегрирования, используемые при решении задач механики; [а^], — матрицы геометрических связей соответственно объемных и поверхностных (и контактных поверхностных) конечных элементов, используемые при решении задач механики; [Ь^г)] — матрицы геометрических связей, позволяющие

учесть приложение внешних сил Я(г) = |Д(г)} в дискретных точках (кя — число дискретных точек); [В(е)] — матрица градиентов конечного элемента для задач механики; кдг — число поверхностных конечных элементов; Ь^ — длина (в общем случае криволинейного) ребра, характеризующая поверхностный конечный элемент (д), вдоль которого задана распределенная нагрузка р = [р}.

Основные процедуры альтернирующего метода Шварца. Как было отмечено выше, альтернирующий метод Шварца является итерационным методом [5-7]. Его суть в рамках конечно-элементной технологии состоит в следующем. Рассмотрим два контактирующих тела А и В. На четных итерациях выполняется коррекция компонент векторов перемещений контактных узлов [ик}(а) и [ик}(^) конечно-элементных моделей тел А и В. Для тела А корректирующие выражения имеют вид

[ик }

_ [^ }?А),т, П = 0;

к}(А),т = \ (тт \2n-1, а2га-1 ( (тт \2n-1 (тт \2n-1 \ „ _ 1 о ( ) }(А)+ а(А)(^к}(Б)>в-[^к}(А),ш) , П = 1 2,

где а2А-^ — итерационный параметр; т (1 ^ т ^ МА) — номер текущего узла, лежащего на контактной поверхности БАА тела А,

здесь МА — число контактных узлов на поверхности [и^б— — вектор перемещения сходственной точки в, лежащей на контактной поверхности тела В. Здесь для простоты принято, что 6га = 0 (см. (13)).

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

На нечетных итерациях проводится коррекция компонент векторов узловых сил и возникающих в контактных узлах конечно-элементных моделей тел А и В. Для тела А корректирующее выражение имеет вид

[^к }2А)~Х = [^к }2А),т - а(А),т([Дк + № }(А),т) , П = 0 1, 12, —

Здесь а(АА) т — итерационный параметр; т (1 ^ т ^ МА) — номер текущего узла, лежащего на контактной поверхности тела А, [Як) 8 — вектор контактной узловой силы в сходственной точке в, лежащей на контактной поверхности тела В. Аналогичное соотношение используется для коррекции компонент векторов узловых сил, возникающих в контактных узлах конечно-элементной модели тела В.

Векторы [Як}(А) и [Як}(£) используются при формировании глобальных векторов узловой нагрузки [Д}(А) и [Я}(в) тел А и В. Затем проводится решение двух подобных задач теории упругости. Таким

образом, данный алгоритм состоит в реализации итерационного процесса поочередного задания на поверхностях контакта и векторов перемещений и и векторов контактных сил и , а также в их соответствующей коррекции с тем, чтобы были выполнены либо силовые контактные условия, если в зоне контакта заданы перемещения, либо кинематические, если заданы контактные силы. Вопросы сходимости подобного типа алгоритмов рассмотрены в работах [3, 4].

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

Пример 1. Рассмотрены два одинаковых однородных изотропных цилиндра, расположенных вертикально и контактирующих по одной из торцевых поверхностей (рис. 1). Нижний торец цилиндра В закреплен в направлении вертикали, на верхний торец цилиндра А действует равномерное давление р =10 МПа. В расчетах использованы конечные элементы первого порядка. Особенностью данной задачи является то, что заранее известно значение контактного давления, которое по величине должно совпадать с давлением, приложенным к цилиндру А. Ниже показаны итерационные приближения контактного давления, иллюстрирующие характер Рис. 1. Схема контактного взаи-сходимости: модействия цилиндров

Число итераций.........

Контактное давление • • • • Относительная ошибка • •

10 15 20 25 30

8,980 9,882 9,970 9,996 9,999 1,02 10-1 1,18-10-2 3 • 10-3 4 • 10-4 1 - 10-4

Пример 2. Рассмотрены два контактирующих цилиндра с разными радиусами (см. рис. 1). Нижний торец цилиндра В закреплен в направлении оси О г, на верхний торец цилиндра А действует равномерное давление р =10 МПа. Использованы конечные элементы второго порядка. На рис. 2 показано поле вертикальных перемещений , на рис. 3 — распределение значений компоненты ог тензора напряжений. Из рисунков следует, что нормальные перемещения и нормаль-

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

Рис. 2. Поле вертикальных перемещений Рис. 3. Поле значений компоненты аг тензора напряжений

Пример 3. Исследовано контактное взаимодействие группы из пяти тел: четыре полых цилиндра с фасками на внешних поверхностях установлены один на другой внутри толстостенной трубы с зазором Ьг = 0,01 мм. Во всех внутренних цилиндрах задано равномерное тепловыделение, которое линейно изменялось во времени до предельного значения ду = 45 кВт/м. Температура наружной поверхности внешнего цилиндра (толстостенной трубы) полагалась постоянной: = 623 ^ Нижние торцы внешней трубы и нижнего внутреннего цилиндра закреплены в направлении оси Ох. На наружной поверхности внешней трубы задано давление рг = 16 МПа. На верхней торцевой поверхности верхнего цилиндра задано давление рг = 2 МПа. Постановка задачи рассмотрена в работе [13].

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

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

Рис. 4. Температурное поле контактирующих тел во всей области (а) и для двух верхних цилиндров (б)

Рис. 5. Поля радиальных (а) и осевых (б) перемещений

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

Рис. 6. Поля значений компонент радиальных (а) и осевых (б) напряжений

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

Работа выполнена при частичной финансовой поддержке гранта Президента РФ для государственной поддержки ведущих научных школ (проект НШ-255.2012.8) и поддержке топливной компании "ТВЭЛ" (Госкорпорация "Росатом").

ЛИТЕРАТУРА

[1] Сакало В.И., Косов В.С. Контактные задачи железнодорожного транспорта. Москва, Машиностроение, 2004, 496 с.

[2] Бабин А.П., Зернин М.В. Конечно-элементное моделирование контактного взаимодействия с использованием положений механики контактной псевдосреды. Известия РАН. Механика твердого тела, 2009, № 4, с. 84-107.

[3] Цвик Л.Б. Принцип поочередной непрерывности при решении задач теории поля по частям. Докл. АН СССР, 1978, т. 243, вып. 1, с. 74-77.

[4] Цвик Л.Б. Принцип поочередности в задачах о сопряжении и контакте твердых деформируемых тел. Прикладная механика, 1980, т. 16, № 1, с. 13-18.

[5] Станкевич И.В., Яковлев М.Е., Си Ту Хтет. Разработка алгоритма контактного взаимодействия на основе альтернирующего метода Шварца. Вестник МГТУ им. Н.Э. Баумана. Сер. Естественные науки, 2011, спец. вып. Прикладная математика, с. 134-141.

[6] Станкевич И.В., Яковлев М.Е., Си Ту Хтет. Математическое моделирование контактного взаимодействия упругопластических сред. Наука и образование: электронное научно-техническое издание, 2012, № 4. URL: http://technomag.edu.ru/doc/353180.html (дата обращения 04.04.2012).

[7] Яковлев М.Е. Математическое моделирование поликонтактного взаимодействия. Вестник МГТУ им. Н.Э. Баумана. Сер. Естественные науки, 2012, спец. вып. № 2 Математическое моделирование в технике, с. 219-224.

[8] Шлыков Ю.П., Ганин Е.А., Царевский С.Н. Контактное термическое сопротивление. Москва, Энергия, 1977, 328 с.

[9] Котович А.В., Станкевич И.В. Решение задач теории упругости методом конечных элементов. Москва, Изд-во МГТУ им. Н.Э. Баумана, 2012, 106 с.

[10] Зарубин В.С., Станкевич И.В. Расчет теплонапряженных конструкций. Москва, Машиностроение, 2005, 352 с.

[11] Галанин М.П., Савенков Е.Б. Методы численного анализа математических моделей. Москва, Изд-во МГТУ им. Н.Э. Баумана, 2010, 591 с.

[12] Котович А.В., Станкевич И.В. Решение задач теплопроводности методом конечных элементов. Москва, Изд-во МГТУ им. Н.Э. Баумана, 2010, 84 с.

[13] Богатырь С.М., Галанин М.П., Горбунов-Посадов М.М., Гусев А.С., Еременко А.С., Ермаков А.В., Кузнецов В.И., Лукин В.В., Новиков В.В., Родин А.С., Салатов А.В., Сыпченко М.В., Фальков А.А., Шаповалов К.Л. Комплекс программ для вероятностных расчетов термомеханики тепловыделяющих элементов. Вестник МГТУ им. Н.Э. Баумана. Сер. Естественные науки, 2012, спец. вып. № 1 Прикладная математика и механика, с. 66-75.

Статья поступила в редакцию 15.05.2013

Ссылку на эту статью просим оформлять следующим образом:

Богатырь С.М., Галанин М.П., Кузнецов В.И., Новиков В.В., Родин А.С., Яковлев М.Е., Крупкин А.В., Лукин В.В., Станкевич И.В. Математическое моделирование термоупругого контактного взаимодействия осесимметричных тел. Инженерный журнал: наука и инновации, 2013, вып. 4. URL: http://engjournal.ru/catalog/mathmodel/ hidden/667.html

Богатырь Сергей Михайлович — нач. отдела ОАО «ВНИИНМ им. акад. А.А. Бочвара».

Галанин Михаил Павлович — д-р физ.-мат. наук, зав. отделом Института прикладной математики им. М.В. Келдыша РАН, проф. кафедры «Прикладная математика» МГТУ им. Н.Э. Баумана. e-mail: galan@keldysh.ru

Кузнецов Владимир Иванович — нач. лаборатории ОАО «ВНИИНМ им. акад. А.А. Бочвара».

Новиков Владимир Владимирович — канд. техн. наук, зам. генерального директора ОАО «ВНИИНМ им. акад. А.А. Бочвара».

Родин Александр Сергеевич — канд. физ.-мат. наук, науч. сотруд. Института прикладной математики им. М.В. Келдыша РАН.

Яковлев Максим Евгеньевич — аспирант, асс. кафедры «Прикладная математика» МГТУ им. Н.Э. Баумана. e-mail: me-yakovlev@rambler.ru

Крупкин Антон Владимирович — ст. науч. сотруд. ОАО «ВНИИНМ им. акад. А.А. Бочвара».

Лукин Владимир Владимирович — канд. физ.-мат. наук, науч. сотруд. Института прикладной математики им. М.В. Келдыша РАН, доц. кафедры "Прикладная математика" МГТУ им. Н.Э. Баумана. e-mail: vvlukin@gmail.com

Станкевич Игорь Васильевич — д-р техн. наук, проф. кафедры «Прикладная математика» МГТУ им. Н.Э. Баумана. e-mail: aplmex@yandex.ru