Численное решение контактных задач с учетом деформации ползучести Текст научной статьи по специальности «Физика»

УДК 539.3

И. В. Станкевич

ЧИСЛЕННОЕ РЕШЕНИЕ КОНТАКТНЫХ ЗАДАЧ С УЧЕТОМ ДЕФОРМАЦИИ ПОЛЗУЧЕСТИ

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

E-mail: aplmex@yandex.ru

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

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

Математическая постановка контактной задачи МДТТ. Рассмотрим два трехмерных однородных и изотропных вязкоупругих контактирующих тела A и B, занимающих в пространстве R3 области Ga и GB и ограниченных кусочно-гладкими границами OGa и dGB. Введем декартовую систему координат OXiX2X3 (см. рисунок).

Математическая формулировка квазистатической несвязанной контактной задачи МДТТ в этом случае будет включать: • уравнения равновесия

j{u,T} + Qi(X,t) = 0, X G Ga; i,j = 173; a G{A,B}; (1)

Схема контактного взаимодействия двух тел

• граничные условия (кинематические и силовые)

и(Х,*)|5? = и0(X,*), X е Б? С два; а е{А,В}; (2) ац{и,Т}пц\3а = р(X,*), X е Б? С два;

1,3 = ТД а е{А,В}; (3)

• соотношения Коши (безындексная запись)

£Ф,*) = Ле/u(X,*); X е ; а е{А,В}; (4)

• определяющие уравнения (безындексная запись)

а = Е (е,Т); (5)

где и^, *) — вектор перемещения точки, определяемой радиусом-вектором X = хе^; и0(X,*) — вектор перемещения точки, расположенной на поверхности Б®, а е {А, В}; £ = £це^ 0ец — тензор деформации; а = аце^ 0 ец — тензор напряжений; Q(X,t) = Qi (X, —

вектор массовых сил, здесь X е а е {А, В}; р^, *) = р^(X, — вектор внешней нагрузки, действующей на поверхности , здесь

V

X е Б® С ; а е {А, В}; Е(£,Т) - тензор-оператор определяющих уравнений.

Так как контактная задача рассматривается в квазистатической постановке, поэтому время * является параметром задачи и изменяется в пределах 0 < * < *тах, где *тах — максимальное (общее) время процесса деформирования.

Запись ац{и, Т} означает, что выполнены все необходимые преобразования, позволяющие тензор напряжений а выразить через век-

тор перемещений и(Х,£) и температуру Т(Х,£). В дальнейшем аргументы у функций будем опускать, если это не будет приводить к противоречию. Кроме массовых Q(X, ¿) и поверхностных р(Х,^) сил в отдельных и°° точках внешней поверхности дСа (а е {А, В}) с фиксированными координатами Хк (к = 1,и°°) могут быть заданы дискретные силы Ик (£).

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

и£(Х,*) - иВ(X, ¿) = 8п(X), X е Бк; (6)

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

аА^Н -оВ(X,*) < 0, X е Бк; (7)

А В

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

т~> А В « А В

В; оА, оВ — проекции векторов напряжений оА и оВ на направление внешней нормали к границе тела А. Здесь Б^ и Б^ и Б^ = дСа, шев(Б? П Б?) = 0, шев(Б? П Б£) = 0, ше8(Б£ П Б£) = 0, а е {А, В}.

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

Для решения задачи (1)-(7) необходимо иметь заданный в явном

V

виде тензор-оператор определяющих уравнений Е(е,Т). Структура

определяющих уравнений (5) зависит от программ нагружения и проявляющихся при этом свойств материала. При решении задач МДТТ с учетом деформации ползучести широкое применение нашли различные варианты теории наследственной ползучести и трех основных технических теорий - старения, течения и упрочнения [2-5]. Известны также теории, использующие для описания ползучести аппарат структурных моделей и механических аналогов [6-8].

Большинство теорий ползучести удовлетворительно описывает деформирование при постоянных и/или медленно изменяющихся нагрузках. Анализ напряженно-деформированного состояния при переменных нагрузках лучше описывается с использованием теорий течения и упрочнения, при этом теория упрочнения имеет некоторые преимущества перед теорией течения, так как несколько точнее аппроксимирует результаты экспериментов [2]. С точки зрения организации вычислительного процесса технические теории имеют известные преимущества перед наследственными [2, 3].

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

Для коррекции кинематических и силовых условий в зоне контакта можно использовать алгоритм, основанный на альтернирующем методе Шварца [9-12].

Основные процедуры альтернирующего метода Шварца. Альтернирующий метод Шварца является итерационным методом. Его суть в рамках конечно-элементной технологии состоит в следующем. На четных итерациях выполняется коррекция компонент векторов перемещений контактных узлов {ик }(а) и {ик }(В) конечно-элементных моделей тел А и В (см. рис. 1). Для тела А корректирующие выражения имеют вид

}Й),т, п = 0;

{и}й),т = | {^^+«га ({и^ - №), (8)

п = Т, 2,...,

где а(А)т — итерационный параметр; т (Т ^ т ^ М4) — номер текущего узла, лежащего на контактной поверхности Б44 тела А, здесь М4 — число контактных узлов на поверхности Б4; {и}(П-1 — вектор перемещения сходственной точки 8, лежащей на контактной поверхности Б^ тела В. Здесь для простоты принято, что 5п = 0 (см. (6)).

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

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

}(А) = {Н-к }(А) ,т а(А),т }(В), в + }(А), т) ,

п = 0, Т, 2, ...; (9)

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

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

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

Таким образом, данный алгоритм состоит в реализации итерационного процесса поочередного задания на поверхностях контакта БА

пЯ « А Я А

и ¿Я векторов перемещений иА и иЯ и векторов контактных сил рА и рЯ, а также в их соответствующей коррекции с тем, чтобы были выполнены либо силовые контактные условия, если в зоне контакта заданы перемещения, либо кинематические, если заданы контактные силы. Вопросы сходимости подобного типа алгоритмов рассмотрены в работах [11, 12]. Примеры численного решения контактных задач теории упругости с помощью данного алгоритма приведены в работах [9, 10].

Основные алгоритмы для решения краевых контактных задач МДТТ с учетом деформации ползучести. При решении МКЭ краевых задач МДТТ с учетом деформаций ползучести весьма часто используют схемы Эйлера — явную или неявную. В зависимости от особенностей рассматриваемой задачи алгоритм решения строится либо в соответствии с методом начальных напряжений, либо — методом начальных деформаций [3, 13]. Метод начальных деформаций при решении задач с учетом ползучести используется чаще, поскольку применение метода начальных напряжений для этого класса задач технически значительно сложнее. Ниже рассматриваются явная и неявная схемы Эйлера в сочетании с МКЭ. Обе схемы формулируются в соответствии с методом начальных деформаций. Определяющее соотношение, для примера, было выбрано в форме теории течения.

Явная схема Эйлера. Данная схема реализуется в соответствии со следующим алгоритмом.

1. В начальный момент времени ¿к = ¿0 = 0 принимается, что во всех конечных элементах (е) конечно-элементных моделей тел А и В деформация ползучести отсутствует {ес(е)(£к)}(а) = {ес(е){(¿0}(а) = {0}, а е {А, В}. В соответствии с процедурами МКЭ формируются глобальные матрицы жесткости [К](а) и глобальные векторы нагрузки {Л}(„), а е {А, В} [1].

2. Решается контактная задача теории упругости [9, 10].

3. Для каждого конечного элемента (е) тел А и В вычисляется вектор напряжений {о(е)}(а), а е {А, В}.

4. Из принятых определяющих соотношений находится вектор скорости деформации ползучести {¿с(е)(4)}(а) = {/({^(е)}(а), )}, а € {А, В}.

5. Для момента времени ¿к+1 вычисляется вектор деформации ползучести

{£с(е)(^+1)}(а) = {ес(£)(4)}(«) + Д*{ёс(е)(4)}(а),а € {А, В}.

где Д£ = ¿к+1 — — шаг по времени.

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

{Д ^О}^) = / [В(е)]т[Н(е)]{£с(е)(^+1)}и^, а € {А, В}.

.¡V (=)

В данном случае фиктивная начальная деформация принимается равной деформации ползучести.

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

{Д£С}(в) = ^[а(е)]т{ЛЙ)(¿к+1)}(а), а € {А, В}.

е=1

8. Формируются суммарные глобальные векторы нагрузки {Де}(а) = {Д}(«) + {Дес}(а), а € {А, В}, учитывающие действие заданной нагрузки ({Д}(а)) и нагрузки, определяемой деформацией ползучести ({Дес}(а)).

9. Решается контактная задача теории упругости.

10. Если заданный интервал времени пройден (4+1 ^ £тах), то вычисления заканчиваются, в противном случае делается переход к пункту 3.

Глобальные матрицы жесткости [К](а) и векторы нагрузки {Д}(а), а € {А, В} строятся на каждом шаге по времени только в том случае, если их компоненты изменяются, например, зависят от температуры, которая, в свою очередь, зависит от времени.

Очевидным преимуществом этой схемы является ее простота и экономичность. Однако явная схема весьма чувствительна к выбору величины шага Д£ по времени.

Неявная схема Эйлера. Данная схема может быть реализована в соответствии со следующим алгоритмом.

1. В начальный момент времени = ¿0 = 0 принимается, что во всех конечных элементах (е) конечно-элементных моделей тел А и В деформация ползучести отсутствует {ес(е)(£к)}(а) = {ес(е) (£0)}(а) = {0},

а е {А, В}. В соответствии с процедурами МКЭ формируются глобальные матрицы жесткости [К](а) и глобальные векторы нагрузки {Д}и, а е {А, В}.

2. Решается контактная задача теории упругости.

3. Для каждого конечного элемента (е) тел А и В вычисляется вектор напряжений {о(е)}(а), а е {А, В}.

4. Из принятых определяющих соотношений находится вектор скорости деформации ползучести {¿с(е)(4)}(а) = {/({а(е)}(а),¿к)}, а е {А, В}.

5. Фиксируется вектор {е0и(*к+1)}(а) = {ес(е)(4 )}(а) и принимается, что {¿с(е) (*к+1)}(а) = {£Гс(е)(^к)}(а), а е {А, В}.

6. Для момента времени ¿к+1 вычисляется вектор деформации ползучести

К(е)(4+1)}(а) =

= К(е)(4)}(а) + Д* [(1 - в^(¿к)}(а) + в{¿Ф)(*к+1)}(а)] ,

где Д* — шаг по времени (Д* = *к+1 — *к); в — параметр (в ^ 0,5), а е {А, В}.

7. Оценивается сходимость итераций внутри рассматриваемого шага по времени: если для тел А и В выполняется неравенство

max

где 8£с — заданная точность, то выполняется переход к пункту 13, в противном случае — к пункту 8.

8. Формируются векторы {Лес}и = Х^=1[а(е)]т{ДЙ)(*к+1)}(а), а е {А, В}, здесь {^(4+1)} (а) = /ум [В(е)]т[Н(е)]{£с(е)(*

9. Решается контактная задача теории упругости.

10. Для каждого конечного элемента (е) тел А и В вычисляется вектор напряжений {о(е)}(а), а е {А, В}.

11. Фиксируются векторы {^(¿к+О}^) = {ес(е)(*к+1)}(а) и вычисляется векторы скоростей деформации ползучести с учетом полученных векторов напряжений {о(е)}(а), (а е {А, В}

{¿с(е)(*к+1)}(а) = {/({о(е)}(«),*к+1)} , а е {А, В}.

12. Выполняется переход к пункту 6.

13. Если заданный интервал времени пройден (¿к+1 ^ ¿таж), то вычисления заканчиваются, в противном случае делается переход к пункту 5.

При реализации неявной схемы глобальные матрицы жесткости [К](а) и векторы нагрузки {Д}(а), а е {А, В}, как и в случае явной схемы, могут быть построены в начале расчета и затем оставаться

постоянными. Если их компоненты изменяются, например, зависят от температуры, которая в свою очередь зависит от времени, то матрицы жесткости [К](а) и векторы нагрузки {Д}(а), а € {А, В}, строятся на каждом шаге по времени, однако на итерациях внутри шагов остаются постоянными.

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

Одной из модификаций как явной, так и неявной схем, направленной на повышение вычислительной эффективности, является следующая [14]. После того, как определены векторы напряжений {^(е)}(«), (а € {А, В}) , вычисляются глобальные векторы невязки узловых сил конечно-элементных моделей тел А и В

{£}а = {Д}а — V[а(е)]т / [В(е)]т{а(е)}аА,, 1=1 •/v(в)

которые затем прибавляются к правым частям основных матричных уравнений [1]

[К](а){и }а = {Д}а + {Д£е }а + 7 {£ }а,

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

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

Работа выполнена в рамках гранта поддержки ведущих научных школ № НШ-255.2012.8.

СПИСОК ЛИТЕРАТУРЫ

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

2. М а л и н и н Н. Н. Прикладная теория пластичности и ползучести. М.: Машиностроение, 1975. 398 с.

3. Б о й л Дж., С п е н с Дж. Анализ напряжений в конструкциях при ползучести. М.: Мир, 1986. 360 с.

4. Ползучесть элементов машиностроительных конструкций / А.Н. Подгорный, В.В. Бортовой, П.П. Гонтаровский и др. Киев: Наукова думка, 1984. 264 с.

5. Ильюшин А. А., Победря Б. Е. Основы математической теории термо-вязкоупругости. М.: Наука, 1970. 280 с.

6. Новожилов В. В., Кадашевич Ю. И. Микронапряжения в конструкционных материалах. Л.: Машиностроение, 1990. 223 с.

7. Г о х ф е л ь д Д. А., С а д а ко в О. С. Пластичность и ползучесть элементов конструкций при повторных нагружениях. М.: Машиностроение, 1984. 256 с.

8. З а р у б и н В. С. Прикладные задачи термопрочности элементов конструкций. М.: Машиностроение, 1985. 296 с.

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

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

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

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

13. Ц в и к Л. Б. Принцип поочередности в задачах о сопряжении и контакте твердых деформируемых тел // Прикладная механика. 1980. Т. 16. № 1. С. 1318.

14. Т е м и с Ю. М. Прикладные задачи термопластичности и термоползучести // Машиностроение. Энциклопедия: В 3 т. / Под общ. ред. К.С. Колесникова. М.: Машиностроение, 1994. Т. 2. С. 226-272.

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