CC BY
15
МАТЕМАТИКА И МЕХАНИКА
УДК 539.3
Н.В. Овчинникова, Ю.В. Чеботаревский
О ПРИМЕНЕНИИ УПРОЩЕННОЙ МОДЕЛИ С ЭКВИВАЛЕНТНОЙ НАГРУЗКОЙ
К РЕШЕНИЮ КОНТАКТНОЙ ЗАДАЧИ
Рассматривается контактная задача о напряженно-деформированном состоянии упругопластической среды на основе расчетной схемы, предполагающей замену реального контактного воздействия на материал среды действием эквивалентной ему распределенной нагрузки. Определены рациональные пределы применимости упрощенной модели для решения упругопластических задач в статической и динамической постановках.
Контактная задача, распределенная нагрузка, напряженное состояние, пределы применимости
N.V. Ovchinnikova, Yu.V. Chebotarevsky
SIMPLIFIED EQUIVALENT MODELLING FOR THE ANALYSIS OF CONTACT PROBLEMS
The paper analyzes the contact problem relating the strain-stress state in the elastoplastic domain. The problem can be solved by substitution of the actual contact forces acting on the material medium for the equivalent distributed loading applied to the contact surface. Applicability limits of the simplified model simulated to solve the stress-strain problems under static or dynamic conditions have been determined.
Contact problem, distributed load, deflected mode, limits of applicability
Сложности, возникающие при непосредственном решении контактных задач, привели к появлению приближенных моделей взаимодействия контактирующих тел. Один из таких подходов предложен в [1-3], где решение исходной контактной задачи сводится к исследованию напряженно-деформированного состояния упругопластической среды под действием эквивалентной нагрузки. Эффективность ее применения к решению статических задач показана в [1-3]. Однако целесообразность ее применения для решения упругих и упругопластических задач в динамической постановке не изучена.
Целью настоящей работы является определение рациональных пределов применимости упрощенной модели для решения упругопластических задач в статической и динамической постановках с точки зрения уменьшения затрат машинного времени на проведение вычислений и снижения технических требований к используемой вычислительной аппаратуре по сравнению с затратами на получение аналогичных решений контактных задач.
Учитывая, что воздействие индентора согласно рассматриваемой модели носит локальный характер, выделим в полубесконечной среде некоторый конечный объем радиусом R и толщиной h. Его размеры должны быть такими, чтобы за пределами выделенного объема среда находилась в невозмущенном состоянии. Тогда на каждом этапе вычислений исследование осесимметричного напряженно-деформированного состояния среды под действием эквивалентной нагрузки, отображающей воздействие индентора, в цилиндрической системе координат может быть сведено аналогично тому, как это было сделано в [4-6], к отысканию поля скоростей точек среды, обеспечивающего безусловный минимум функционала Ф(г, z, t):
5Ф(г, г) = 11 ( огг (Буг )+ 5уг ^ + с* (Ч)+ 5уг р ^+
э э ду Л а (1)
+ с22—^ )+ сгг — (5у2 )+ 5у2р —^ rdzdr - } 5у24(г, г)гШг = 0, Эг Эг Эг ) 0
где Су (/, j = г, 0, г) - компоненты тензора напряжений; уг, 5уг, - компоненты вектора скорости точек среды и их вариации; а и д(г, г) - соответственно радиус пятна контакта и эквивалентная нагрузка, определяемая путем решения задачи о движении индентора; р - плотность материала среды.
Разобьем выделенный объем среды на элементы, связанные между собой в п узловых точках, и выразим составляющие вектора скорости движения произвольной точки среды у,(г, г, г) (/ = г, г) через составляющие вектора скорости узловых точек ук(г) следующим образом:
п
V (г, 2, г )= X N (Г, Ук1 (г ) = Ык (Г, Ук1 (г) (/ = г г). (2)
к=1
Подставляя (2) в уравнение (1), после выполнения ряда преобразований, аналогичных проделанным в [4-6], получаем
п [ Кк ( Э / \ Д/\П
5П = X5^Ш с^К)+ N.+ )+рX^каУтг
гш2шг} +
(3)
п I к2 "2 / ^ / \ ^ I \ п 1 а |
+ ^ ^ 1 /( ^ ^к К ^к )+Р Х "к^ ^ - 0 ^ * = 0.
к=1 [00 V Эг г эг т=1 dг
п [К2 к2 ( эил э1\ Л а
-1 \ к / ' ~ гг \ ' к I ' г ¿—I ~ ' к~ ' т ,
к=1 [ 0 0 V Эг Эг т=1 dг
Отсюда в силу независимости и произвольности вариаций скоростей 5укг и приходим к следующей системе обыкновенных дифференциальных уравнений первого порядка для определения скоростей в узловых точках разбиения среды на элементы:
X = - К12 1 ( с гг N) + N. ^ + с гг А N )) г dzdr;
т=1 dг 0 0 V Эг г Эг ')
= -Кк(о»)+^г"Э"^)! Гdzdг + а^д(г,г,г)гdг , (4)
т=1 ш пп V Эг Эг ) 0
т=1 00\ Ш иг у 0
где элементы матрицы масс определяются соотношением
Кк
N А
1 т
Мкт =Р 11 NkNmrdrdz . (5)
00
При проведении конкретных вычислений уравнения (4) должны быть дополнены физическими уравнениями, соответствующими характеру поведения материала среды в процессах нагружения и разгрузки, и геометрическими соотношениями Коши аналогично тому, как это сделано в [4-6].
Сформулированная выше конечно-элементная схема численного решения поставленной задачи реализована с использованием программного комплекса ABAQUS. Расчеты проводились с применением треугольных конечных элементов и квадратичных функций формы на основе метода Ньютона.
Прежде чем заняться определением этих границ применимости рассматриваемой модели с точки зрения поставленной цели, следует убедиться в том, что построенная конечно-элементная схема действительно отображает исследуемый процесс и получаемые на её основе результаты достоверны. Этот вопрос здесь решается двумя способами:
- путём сравнения результатов, полученных на базе разработанной модели, с результатами известного аналитического решения А.Н. Динника, полученного для случая статического нагружения двух соприкасающихся тел в виде полусферы и полупространства с плоской ограничивающей его поверхностью;
- путем их сравнения с результатами численных решений, полученных по методике, разработанной в [4-6], для различных вариантов нагружения индентора в случае классической постановки контактной задачи.
Второй способ позволяет не только оценить достоверность получаемых результатов, но и провести сравнительный анализ затрат вычислительных ресурсов на их получение с применением упомянутой выше методики на основе упрощенной модели.
Для того чтобы задать интенсивность эквивалентной нагрузки, предварительно необходимо решить задачу о движении индентора [1-3]. При упругом поведении материала эта задача может быть рассмотрена отдельно от задачи о напряженно-деформированном состоянии среды. При упругопла-стическом состоянии эти задачи оказываются взаимосвязанными. Для того чтобы облегчить задание эквивалентной нагрузки в процессе решения конечно-элементной задачи, аппроксимируем найденные ее значения рядом Тейлора [1-3]:
Ж ) = Чо 1 + Ё
{ (к)/Л\ л гк
"д(к )(0)
к=1 к! у
(6)
Такой подход позволяет задать эквивалентную нагрузку в аналитической форме и тем самым избежать необходимости формирования массива её значений в отдельных точках контактной поверхности в процессе вычислений.
Хорошее совпадение численных результатов, полученных на основе упрощенной модели, для статического нагружения индентора в пределах, не вызывающих появления в среде пластических деформаций, с результатами аналитического решения [7], вообще говоря, предсказуемо. Так как в основу упрощенной модели при определении равнодействующей сил контактного взаимодействия и эквивалентной нагрузки заложены зависимости между радиусом пятна контакта, распределением давления на контактной поверхности и перемещением центральной точки контакта, полученные на базе упомянутого выше решения.
При статическом нагружении индентора результаты расчета напряженно-деформированного состояния среды как на базе предлагаемой модели, так и на основе аналитического решения, определяются величиной усилия, прикладываемого к индентору, и не зависят от его массы. При динамическом нагружении индентора в пределах упругого поведения материала среды значение его приведенной массы заметно влияет на частоту колебаний точек среды. Поэтому в задаче о движении инденто-ра при подсчете приведенной массы, кроме непосредственно массы индентора, следует учитывать и пришедшую в движение массу среды.
Ниже показаны графики поведения интенсивности напряжений по времени при импульсном нагружении индентора силой Г = Го = 19,62 Н за промежуток времени М = 103 с в точке её наибольших значений. На рисунке б приведены графики, построенные по расчетным данным, полученным на базе предлагаемой модели. На рисунке а - по результатам решения аналогичной контактной задачи способом, разработанным в [4-6]. Расчеты проведены для значений параметров (Е = 2,1-Ю11 Па, V = 0,3).
а б
Графики поведения интенсивности напряжений с течением времени
Сравнительный анализ приведенных графиков показывает их хорошее совпадение. Погрешность вычислений при использовании упрощенной модели для наибольших значений интенсивности напряжений не превышает 5 %. Выигрыш во времени обусловлен возможностью решения задачи о движении индентора независимо от движения среды. Затраты машинного времени ниже не менее чем на 5 % по сравнению со временем, необходимым для получения решения контактной задачи в полном объеме. Кроме того, для её численной реализации возможно применение компьютерной техники со значительно меньшим (на порядок) вычислительным ресурсом.
При появлении в среде пластических деформаций задача о движении индентора не может рассматриваться отдельно от задачи о движении среды. Обе задачи должны решаться совместно ввиду необходимости реализации итерационного процесса по определению секущего модуля упругости в задаче о движении индентора по результатам расчета напряженно-деформированного состояния среды на каждом шаге вычислений по формулам:
М I \ 1 П а(т"!) , ,
£о) = £, Е{гп) = 1 £ ^ (т = 1,2,3,...) . (7)
nk=1 е[т 1
В (7) ^ и $еи - значения интенсивности напряжений и интенсивности полной деформации в
точках возмущенной области среды соответственно, п - их количество, а т - номер итерации. Легко убедиться, что в частном случае при идеально упругом поведении материала среды секущий модуль, определяемый по формуле (7), становится равным модулю упругости её материала.
В табл. 1 в качестве примера приведено сопоставление расчетных данных, полученных на базе предлагаемой модели и способом, разработанным в [4-6], в случае нагружения индентора статическим усилием ¥ = ¥0 = 39,24 Н. Результаты проведенных расчетов показали, что значения интенсивности напряжений, полученных обоими способами, практически совпадают с точностью до второго знака после запятой. Расхождение численных значений интенсивности пластических деформаций в точках их наибольших значений не превышает 1 %. Наибольшее несовпадение расчетных данных интенсивности пластических деформаций имеет место в граничных точках пластической зоны, но оно несущественно в силу малости самих значений.
Таблица 1
Координата 2 точки, ■10-4 м е и 104 а и 10-8 Па
из решения контактной задачи из решения на основе модели из решения контактной задачи из решения на основе модели
0 0 0 0,517377 0,497223
0,125 0 0 1,06074 1,07956
0,25 0 0 1,54985 1,61685
0,375 0,233802 0,322716 1,85379 1,8961
0,5 1,14673 1,53697 2,0628 2,06453
0,625 2,51198 2,74618 2,06498 2,071
0,75 3,6547 3,83884 2,09431 2,09915
0,875 4,0772 4,1924 2,10474 2,10772
0,1 4,48226 4,51488 2,11517 2,11611
1,125 4,39502 4,3607 2,1128 2,11196
1,25 4,2633 4,18534 2,10975 2,10783
1,375 3,77971 3,72015 2,09728 2,0958
1,5 3,30028 3,25571 2,08547 2,08434
1,625 2,67092 2,60147 2,06919 2,06734
1,75 2,00645 1,9456 2,05256 2,05095
1,875 1,26473 1,25969 2,03313 2,03297
2 0,283531 0,385823 2,05523 2,04851
2,125 0,0313116 0,0889457 1,9644 1,97092
Для обеспечения удовлетворительного совпадения получаемых результатов в приведенном выше примере при нагружении усилием ¥ = ¥0 = 39,24 Н, потребовалось проведение восьми итераций
(табл. 2). При дальнейшем увеличении числа итераций отличие последующих значений секущего модуля упругости от предыдущих не превышает 2,2 %. С возрастанием величины усилия совпадение результатов расчета, полученных обоими способами, улучшается, но количество итераций, необходимых для обеспечения заданной точности вычислений, при этом увеличивается.
Таблица 2
Итерация Секущий модуль Е^"^ .-1011 Па ер в точке с координатами г = 0, г = 110-3 м, 104 Си в точке с координатами г = 0, г = 110-3 м, 108 Па
0 2,1 8,0642 2,20415
1 1,6589 2,96187 2,0768
2 1,9413 6,16109 2,15726
3 1,744 3,8692 2,09983
4 1,8816 5,43049 2,13906
5 1,7751 4,19835 2,10814
6 1,8405 4,942 2,12684
7 1,8032 4,51488 2,11611
Из вышесказанного следует, что для расчета напряженно-деформированного состояния упру-гопластической среды под действием эквивалентной статической нагрузки технические требования к вычислительным инструментам те же, что и при решении задачи в полном объеме, а затраты машинного времени в среднем одинаковы.
Очевидно, что с увеличением количества итераций время счета возрастает практически пропорционально их числу. При использовании упрощенной модели в случае приложения к индентору динамической нагрузки реализация итерационного процесса необходима на каждом шаге вычислений по времени. Поэтому по результатам анализа процедуры расчета напряженно деформированного состояния среды на базе упрощенной модели в динамической постановке сделан вывод о нецелесообразности её применения для проведения оперативных расчетов в связи со значительным ростом затрат вычислительных ресурсов из-за необходимости многократного перерасчета секущего модуля в задаче о движении индентора.
ЛИТЕРАТУРА
1. Овчинникова Н.В., Чеботаревский Ю.В. О некоторых особенностях применения метода конечных элементов к решению контактной задачи на базе программного комплекса ABAQUS // Изв. Саратовского университета. Новая серия. Сер. Математика. Механика. Информатика. 2009. Т. 9. Вып. 2. С. 82-88.
2. Овчинникова Н.В., Чеботаревский Ю.В. К расчету напряженно-деформированного состояния упругопластического полупространства, контактирующего с абсолютно жестким индентором // Вестник СГТУ. 2010. № 4 (51). Вып. 3. С. 10-17.
3. Овчинникова Н.В., Чеботаревский Ю.В. О движении абсолютно жесткого индентора, взаимодействующего с упругопластической средой // Известия высших учебных заведений. Поволжский регион. Технические науки. 2010. № 3 (15). С. 151-164.
4. Овчинникова Н.В., Чеботаревский Ю.В. Вариационное уравнение движения континуума «жесткий индентор - деформируемая среда» // Вестник СГТУ. 2011. № 4 (60). Вып. 2. С. 48-57.
5. Овчинникова Н.В., Чеботаревский Ю.В. Применение метода множителей Лагранжа к решению контактной задачи о взаимодействии деформируемой среды с относительно жестким инденто-ром // Вестник СГТУ. 2012. № 4 (68). Вып. 2. С. 36-43.
6. Овчинникова Н.В., Чеботаревский Ю.В. Применение метода штрафных функций к выводу вариационного уравнения движения континуума «индентор - деформируемая среда» // Вестник СамГТУ. Сер. Технические науки. 2013. № 1 (37). С. 127-134.
7. Динник А.Н. Избранные труды. Киев: Изд-во Академии наук Украинской ССР, 1952. Т. I.
152 с.
Овчинникова Наталья Владимировна -
кандидат технических наук, доцент кафедры «Техническая механика и детали машин» Саратовского государственного технического университета имени Гагарина Ю.А.
Чеботаревский Юрий Викторович -
доктор технических наук, профессор кафедры «Прикладная математика и системный анализ» Саратовского государственного технического университета имени Гагарина Ю.А.
Natalya V. Ovchinnikova -
Ph. D., Associate Professor, Department of Technical Mechanics and Machine Elements,
Yuri Gagarin State Technical University of Saratov
Yuri V. Chebotarevsky -
Dr. Sc., Professor
Department of Applied
Mathematics and System Analysis,
Yuri Gagarin State Technical University of Saratov
Статья поступила в редакцию 12.09.15, принята к опубликованию 11.11.15
