CC BY
188
DOI: 10.18698/1812-3368-2016-1-111-121
УДК 539.3:519.6
ВЫБОР ПАРАМЕТРОВ ВЫЧИСЛИТЕЛЬНЫХ АЛГОРИТМОВ ПРИ РЕШЕНИИ ЗАДАЧИ КОНТАКТНОГО ДЕФОРМИРОВАНИЯ ШЕРОХОВАТЫХ ТЕЛ В ANSYS
М.В. Мурашов
МГТУ им. Н.Э. Баумана, Москва, Российская Федерация e-mail: murashov@bmstu.ru
Численное решение задач контактного деформирования двух шероховатых тел затруднено ввиду большого числа параметров вычислительных алгоритмов и особенностей конечно-элементных моделей, значительно влияющих на результаты решения. Упругопластическая деформация материалов и изменение геометрии контактирующих шероховатых поверхностей приводит к нелинейной задаче деформирования, что часто оказывается дополнительным препятствием для успешного проведения расчета. На примере решения контактной задачи для участка шероховатости микронных размеров в ANSYS рассмотрено влияние на определяемую площадь фактического контакта одного из параметров расширенного метода Лагранжа — допуска на проникание. Приведены рекомендации по выбору значений этого параметра, позволяющие избегать физически неправдоподобных решений. Проанализировано влияние шага приложения нагрузки. Рассмотрены различные стратегии изменения шага нагружения в процессе расчета.
Ключевые слова: контактная тепловая проводимость, шероховатость, метод конечных элементов, упругопластическая деформация, ANSYS.
SELECTION OF PARAMETERS FOR COMPUTATIONAL ALGORITHMS WHILE SOLVING THE PROBLEM OF ROUGH BODIES CONTACT WITH ANSYS
M.V. Murashov
Bauman Moscow State Technical University, Moscow, Russian Federation e-mail: murashov@bmstu.ru
The process of finding a numerical solution to contact deformation problems of two rough bodies is difficult due to many parameters of computational algorithms,as well as special features of finite-element models, which significantly affect the results. Both the elastic-plastic deformation of the materials and the geometry change of the contacting rough surfaces result in appearing a nonlinear deformation problem and often cause an extra difficulty in computing. The paper considers the impact of one parameter of the augmented Lagrangian method (a penetration tolerance) on the determined real contact area. The contact problem solution for the rough area of micron sizes with ANSYS is taken as an example. The author gives recommendations for selecting values of this parameter in order to avoid physically improbable solutions. The impact of the load step is analyzed. Different strategies of the load step modifications during the computational process are considered.
Keywords: thermal contact conductance, roughness, finite element method, elastic-plastic deformation, ANSYS.
Решая задачу контактного деформирования шероховатых поверхностей, можно вычислить площадь фактического контакта, которая
является основным неизвестным в трибологии, тепло-, магнито- и электропроводности в составных телах. Особенно это относится к конструкциям приборов точного приборостроения и микроэлектромеханическим устройствам.
Применение метода конечных элементов к решению задач контактного деформирования двух шероховатых тел открывает возможность создания и исследования моделей контакта, учитывающих структуру шероховатости. В частности, мелкие неровности на поверхности элементов шероховатости при наличии размерных эффектов оказывают существенное влияние на процесс контактного взаимодействия [1, 2]. Важнейшим в таком подходе является возможность учитывать изменение свойств материала в процессе контактирования для микробугорков шероховатости [3]. Однако это вызывает значительные трудности, обусловленные нелинейностью задачи вследствие изменения геометрии и упругопластических свойств материалов, а также большим количеством параметров применяемых вычислительных алгоритмов, совокупное влияние которых на результат существенно, но недостаточно изучено. К таким параметрам относятся допуск на проникание, стратегия выбора шага приращения нагрузки, плотность сетки конечных элементов, контактная жесткость, режим функции автоматического регулирования зазора, начальное проникание и др., а также использование контактных элементов высокого порядка, тетраэдральных или гексаэдральных элементов в разбиении, выбор метода определения областей контакта. Ввиду сложности задачи приоритетным можно полагать не исследование взаимовлияния всех параметров, а поиск наиболее удачной их комбинации с позиции определения площади фактического контакта как результата задачи.
Во многих работах по расчету контактов в программном комплексе А^УБ, например, в работе [4], вопросы выбора параметров контактных алгоритмов опускаются или приводятся лишь использованные значения [5]. В ряде работ применительно к контактным задачам расчета макроразмерных деталей приведены рекомендации по выбору допуска на проникание и контактной жесткости, например, для билинейной и мультилинейной моделей упругопластического материала [6, 7]. Однако в указанных работах рассматривался контакт гладких поверхностей и не ставилась цель точно определить область фактического контакта, а находились напряжения и перемещения.
В настоящей статье рассмотрено поведение фактической площади контакта в зависимости от одного из основных параметров контактного алгоритма на базе расширенного метода Лагранжа — коэффициента допуска на проникание. Расчеты проведены с использованием программного комплекса А^УБ. В алгоритме расширенного метода Ла-гранжа в зоне контакта вводятся гипотетические пружинные элементы
с заданной контактной жесткостью, сила реакции одного тела по отношению к другому вычисляется как произведение гипотетического взаимного проникания тел и контактной жесткости. Обычно полагается, что чем ниже значение проникания, тем выше точность результатов, но хуже сходимость итерационного решения. Ограничение на значение проникания выбирает пользователь программного комплекса, который делает этот выбор методом проб и ошибок на основании имеющегося у него опыта. В статье также рассмотрено влияние прибавления нагрузки.
Постановка задачи. Рассмотрим два контактирующих тела N, О объемами Vм, V0 и размерами 22,5 х 22,5 х 15 мкм, ограниченные неподвижными поверхностями 5м и 50. Высота тел 15 мкм выбрана на основе предварительных вычислений так, чтобы на верхней и нижней поверхностях перемещения были распределены достаточно равномерно. Тела имеют гладкие внешние поверхности и неидеальный контакт неровных поверхностей между собой (рис. 1, а). Обозначим через , , , части поверхности = и и и , и 50, 520, , 50 части поверхности 50 = 50 и 520 и 530 и 50. На поверхность тела N действует внешнее давление р = 8 МПа. Поверхность закреплена от перемещений по оси ^ а поверхности и £0 — от перемещений в ортогональных им направлениях. Трением пренебрегаем.
Рис. 1. Геометрическая схема (а) и модель (б) контакта шероховатых тел
Перейдем к индексным обозначениям осей декартовой системы координат, заменив обозначения x, у, z обозначениями £Ь Х2, Хз. Используется теория течения и аддитивный подход к формированию приращений деформаций. Математическая модель представляет собой уравнения равновесия, обобщенный закон Гука, закон течения, соотношения Коши, уравнения совместности полных деформаций, критерий текучести Мизеса, соотношение для расчета контактных давлений расширенного метода Лагранжа, граничные условия:
Tij,j = £j = e E^ijакк' d£Pj = Sij£ij = 2
d2£ik 92£ji д2£ц 02£jk + т.—tT— = ~—~--Ь j
8x3 6x1 дхгдхк дх3 дхк дх^дхг
2 2 2 ( Р \ 2
- ^2) + (^2 - аз) + (аз - а 1) = 2Ф (в ) ;
р (хг) = Кб + А, хг Е Б°;
азз = Р, хз Е Б1 ,
и1 (х2) = 0, и2 (х1) = 0, х1, х2 Е Б? + Б°,
и (хз) = 0, хз Е Б°,
где а^ и вц — декартовы компоненты тензоров напряжений и деформаций; и — компоненты вектора перемещений; E — модуль Юнга; V — коэффициент Пуассона; бц — символ Кронекера; % — девиа-тор текущих напряжений; а1, а2, аз — главные напряжения; ¿А, А — множители Лагранжа; Ф (вР) — задаваемая функция кривой упрочнения материала; p — контактное давление; K — контактная жесткость; б — зазор между контактирующими поверхностями; азз — компонента тензора напряжений.
Контактирующие поверхности сформированы следующим образом. На ровные поверхности тел, параллельные плоскости xy, наносится сетка узлов размером 4x4, которые затем сдвигаются по нормали к поверхности на 0... 5 мкм в соответствии с масштабированной фрактальной функцией Веерштрасса-Мандельброта [8]. Принятые параметры функции: размер поверхности Ь = 0,0225 м; фрактальная шероховатость С = 1 • Ю-6 мм; фрактальная размерность Б = 2,3; масштабный параметр 7 = 1,5; число волн М = 10; максимальный номер частоты птах = 17; длина отсечки Ь3 = 1 • 10-2 мм. Полученная негладкая поверхность, использованная в расчетах, показана на рис. 1, б.
Задача деформирования решается в квазистатической постановке. Материал тел — медь М2 с модулем упругости Е = 120 ГПа и коэффициентом Пуассона ц = 0,38. Предположим, что материал изотропно упрочняющийся, диаграмма пластического деформирования, взятая
48,1108 144,332 240,554 336,775 0 96,2215 192,443 288,665 384,886
Рис. 2. Поле эквивалентных напряжений, МПа, и форма контактирующих тел до (а) и после (б) деформации
из работы [9], аппроксимировалась мультилинейной кривой с пределом текучести <гт = 69,6 МПа. Для моделирования объема материала использовались изопараметрические конечные элементы SOLID227. Элементами CONTA174 покрыта нижняя поверхность, а TARGE170 — верхняя. Форма контактирующих тел до и после деформации, а также рассчитанное поле напряжений после деформации показаны на рис. 2.
Влияние допуска на проникание. Допуск на проникание обозначается в программном комплексе ANSYS как PTOLN, равен произведению коэффициента допуска на проникание (FTOLN) и средней глубины контактирующих объемных конечных элементов типа SOLID187 или SOLID227 (параметр Average contact pair depth). Максимальное допустимое значение PTOLN не должно превышать размер контактирующих элементов SOLID187 (или SOLID227), т.е. допуск FTOLN не должен быть больше 1. Коэффициент допуска на проникание (FTOLN) по умолчанию имеет значение 0,1. Рекомендуется использовать максимальное значение из диапазона 0,1... 1,0, при котором возможен расчет. В случае повышения допуска на проникание при некоторых его значениях невозможно проведение расчета вследствие непреодолимого прерывания по ошибке искажения элементов (One or more elements have become highly distorted). Больший допуск на проникание позволяет узлам глубже проникать в материал, поэтому менее вероятной становится ситуация появления мельчайших зазоров в результате искривления элементов. Поскольку к области непосредственного контакта (обозначаемой sliding) относятся доли конечных элементов, прилежащие к узлам с положительным прониканием, число элементов
непосредственного контакта возрастает с соответствующим увеличением площади фактического контакта.
В случае достаточно подробной сетки конечных элементов контактирующих тел увеличение БТОЬК слабее изменяет площадь фактического контакта. Если разбиение неудачное, то повышение БТОЬК может привести к увеличению площади контакта в несколько раз (рис. 3, рис. 4).
При БТОЬК < 1 на кривой (см. рис. 3) имеют место участки уменьшения площади фактического контакта, т.е. при наращивании нагруз-
0 12345678 Давление, МПа
в
Рис. 3. Зависимость площади фактического контакта от давления:
а — при увеличении БТОЬК и условии разной поверхностной плотности разбиения сетки конечных элементов (1 — 0,1; 2 — 0,1 для мелкой сетки; 3 — 1,0; 4 — 1,0 для мелкой сетки); б — при постепенном увеличении БТОЬК (1 — 0,1; 2 — 0,5; 3 — 0,9; 4 — 1,0); в — при различном числе шагов нагружения (1 — 88 шагов при БТОЬК = 0,1; 2 — 42 шага при БТОЬК=1,0; 3, 4 — 200 шагов для БТОЬК = 0,1 и БТОЬК=1,0)
Рис. 4. Области непосредственно контакта тела О при р = 0,5 МПа (крупная сетка) и FTOLN = 0,1 (а) и FTOLN = 1,0 (б)
ки некоторые элементы выходят из контакта. Такое поведение обусловлено работой используемого нелинейного контактного алгоритма. При FTOLN = 1 участки спада отсутствуют, что является уже физически правдоподобным. Повышение значения FTOLN заметно сокращает требуемое для расчета общее число итераций и кривая зависимости площади фактического контакта от нагрузки становится более гладкой.
При резком различии размеров поверхностных конечных элементов у контактирующих тел использование значений FTOLN, близких к 1, приводит к искривлению элементов (и останову расчета по соответствующим сообщениям об ошибках), так как средняя глубина контактирующих элементов много больше размера имеющихся мелких поверхностных элементов. В таком случае выполнить расчет возможно только при FTOLN < 1. Это не должно заметно уменьшать площадь фактического контакта в области мелких поверхностных конечных элементов, поскольку PTOLN все равно больше их размера. В областях более грубого разбиения на большие конечные элементы уменьшение площади фактического контакта может быть заметным.
Следует также контролировать величину проникания, так как она является составляющей погрешности определения зазора в контакте и, соответственно, передачи теплоты через зазор.
Влияние выбора шага нагружения. При решении данного класса задач шаг нагружения необходимо задавать переменным. Как правило, число шагов выбирается минимальным, так как на точность это заметно не влияет при достаточно мелком разбиении на конечные элементы. Чаще всего сообщения об ошибках возникают на начальных шагах расчета. В программном комплексе ANSYS задается начальное и минимальное число шагов (number of substeps, minimal number of substeps). При расчете с переменным шагом нагружения программный комплекс ANSYS пытается уложиться в число шагов, указанное как минимальное, начиная с дробления нагрузки на число шагов, заданное
как начальное. Для исключения сообщений об ошибках, возникающих в начале расчета, можно увеличивать начальное число шагов, не меняя минимального числа шагов. Если сообщения об ошибках появляются не на начальном этапе расчета, то необходимо увеличивать минимальное число шагов. Минимальное число шагов должно всегда быть меньше начального числа шагов. При расчете на большие давления в моменты появления новых пятен контакта возникают трудности, как и при первом контакте в начале расчета. Иногда в таких случаях можно пробовать запускать задачу с постоянным числом шагов нагружения.
При недостаточном числе шагов нагружения сам процесс расчета может выполняться без сообщений об ошибках, однако площадь фактического контакта на некоторых шагах расчета будет определена с погрешностью до 70% и более. Увеличение числа шагов нагружения позволяет избавиться от колебаний фактической площади контакта при росте нагрузки (см. рис.3,в). Однако существует предельное число шагов нагружения, превышение которого не ведет к изменению в зависимости. Например, при числе шагов 200 и 500 зависимости практически совпадают.
Зависимость площади фактического контакта от давления представляет собой дискретную функцию вследствие того, что при квазистатическом нагружении нагрузка возрастает до максимальной за заданное число шагов. Площадь фактического контакта меняется скачкообразно, и скачок равен площади очередного поверхностного конечного элемента, входящего в контакт.
При достаточно большом числе шагов нагружения (см. рис. 3, в) все еще заметно некоторое колебание площади фактического контакта при возрастании нагрузки. Причина таких колебаний численная, возможно, связанная с контактным алгоритмом и дискретизацией пространства. В этом случае колебания уже не зависят от числа шагов нагружения, но их значение может составлять до 40% фактической площади. Измельчение сетки помогает избавиться от этих колебаний, однако на практике такой подход крайне затруднителен ввиду ограниченности вычислительных ресурсов.
В случае невозможности дальнейшего деления сетки для повышения точности расчета площади фактического контакта предлагается поступать следующим образом. Зависимость площади контакта от нагрузки обычно имеет линейный характер для малых давлений, хотя из теории упругости следует, что она должна быть степенной с показателем степени меньше единицы [10]. Линейная зависимость получена и при моделировании наноконтактов методами молекулярной динамики для удельных давлений р/Е до 0,05 [11]. Обсуждаемая в литературе и в целом общепринятая линейная зависимость прослеживается и
для зависимостей с мелкими шагами нагружения (см. рис. 3, б). Основываясь на этом для малых давлений (для контакта медь М2-медь М2 p < 50 МПа) следует проводить расчет на максимально возможной мелкой сетке в как можно большем интервале давлений, а затем выполнять линейную экстраполяцию полученной кривой.
Заключение. Коэффициент допуска на проникание, как и шаг увеличения нагрузки в контактных расчетах являются параметрами, назначаемыми пользователем. Проведенные расчеты показали, что при неправильном выборе любого из этих параметров расчет либо невозможен вследствие возникающих сообщений об ошибках, либо вычисленная площадь фактического контакта может иметь погрешность более 70%. Анализ результатов численных экспериментов позволил сформулировать рекомендации, позволяющие проводить расчеты контакта шероховатых тел и существенно повысить точность вычисления площади фактического контакта.
ЛИТЕРАТУРА
1. Murashov M.V., Panin S.D. Modeling of thermal contact conductance // Proceedings of the International heat transfer conference IHTC-14. August 8-13. 2010. Washington, DC, USA. Vol. 6. P. 387-392. DOI: 10.1115/IHTC14-22616
2. Мурашов М.В., Панин С.Д. Моделирование термического контактного сопротивления // Труды пятой Российской национальной конференции по теплообмену. Москва, 25-29 октября 2010. М.: Издательский дом МЭИ. 2010. Т. 7. С. 142-145.
3. Murashov M.V, Panin S.D. Numerical modelling of contact heat transfer problem with work hardened rough surfaces // International Journal of Heat and Mass Transfer. 2015. Vol. 90. P. 72-80. DOI: 10.1016/j.ijheatmasstransfer.2015.06.024
4. Cheng L., Hannaford B. Finite Element Analysis for evaluating liver tissue damage due to mechanical compression//Journal of Biomechanics. 2015. Vol. 48. P. 948-955. DOI: 10.1016/j.jbiomech.2015.02.014
5. Alvarez-Arenal A., Lasheras F.S., Fernandez E.M., Gonzalez I. A jaw model for the study of the mandibular flexure taking into account the anisotropy of the bone // Mathematical and Computer Modelling. 2009. Vol. 50. P. 695-704.
DOI: 10.1016/j.mcm.2008.12.029
6. Lanoue F., Vadean A., Sanschagrin B. Finite element analysis and contact modelling considerations of interference fits for fretting fatigue strength calculations // Simulation Modelling Practice and Theory. 2009. Vol. 17. P. 1587-1602.
DOI: 10.1016/j.simpat.2009.06.017
7. Diaz J.J. del C., Nieto P.J.G., Biempica C.B., Rougeot G.F. Non-linear analysis of unbolted base plates by the FEM and experimental validation // Thin-Walled Structures. 2006. Vol. 44. P. 529-541. DOI: 10.1016/j.tws.2006.04.008
8. Yan W., Komvopoulos K. Contact analysis of elastic-plastic fractal surfaces // Journal of Applied Physics. 1998. Vol. 84. P. 3617-3624. DOI: 10.1063/1.368536
9. Сопротивление деформации и пластичность металлов при обработке давлением / Ю.Г. Калпин, В.И. Перфилов, П.А. Петров, В.А. Рябов, Ю.К. Филиппов. М.: Машиностроение, 2011. 244 с.
10. Шлыков Ю.П., Ганин Е.А., Царевский С.Н. Контактное термическое сопротивление. М.: Энергия, 1977. 328 с.
11. Yang C., Persson B.N.J.Contact mechanics: contact area and interfacial separation from small contact to full contact // Journal of Physics: Condensed Matter. 2008. Vol. 20. No. 21. 215214. 13 p. DOI:10.1088/0953-8984/20/21/215214
REFERENCES
[1] Murashov M.V., Panin S.D. Modeling of thermal contact conductance. Proc. of the International heat transfer conference IHTC-14. August 8-13, 2010, Washington, DC, USA, vol. 6, pp. 387-392. DOI: 10.1115/IHTC14-22616
[2] Murashov M.V., Panin S.D. Modelling of thermal contact resistance. Tr. 5-y Ross. nacionalnoy konf. po teploobmenu [Proc. of the fifth Russian national conf. on heat transfer]. Moscow, MEI Publ., 2010, vol. 7, pp. 142-145 (in Russ.).
[3] Murashov M.V., Panin S.D. Numerical modelling of contact heat transfer problem with work hardened rough surfaces. International J. of Heat and Mass Transfer, 2015, vol. 90, pp. 72-80. DOI: 10.1016/j.ijheatmasstransfer.2015.06.024
[4] Cheng L., Hannaford B. Finite Element Analysis for evaluating liver tissue damage due to mechanical compression. J. of Biomechanics, 2015, vol. 48, pp. 948-955. DOI:10.1016/j.jbiomech.2015.02.014
[5] Alvarez-Arenal A., Lasheras F.S., Fernandez E.M., Gonzalez I. A jaw model for the study of the mandibular flexure taking into account the anisotropy of the bone. Mathematical and Computer Modelling, 2009, vol. 50, pp. 695-704. DOI:10.1016/j.mcm.2008.12.029
[6] Lanoue F., Vadean A., Sanschagrin B. Finite element analysis and contact modelling considerations of interference fits for fretting fatigue strength calculations. Simulation Modelling Practice and Theory, 2009, vol. 17, pp. 1587-1602. DOI:10.1016/j.simpat.2009.06.017
[7] Diaz J.J. del C., Nieto P.J.G., Biempica C.B., Rougeot G.F. Non-linear analysis of unbolted base plates by the FEM and experimental validation. Thin-Walled Structures, 2006, vol. 44, pp. 529-541. DOI:10.1016/j.tws.2006.04.008
[8] Yan W., Komvopoulos K. Contact analysis of elastic-plastic fractal surfaces. J. of Applied Physics, 1998, vol. 84, pp. 3617-3624. DOI: 10.1063/1.368536
[9] Kalpin Yu.G., Perfilov V.I., Petrov P.A., Ryabov V.A., Filippov Yu.K. Soprotivlenie deformatsii i plastichnost metallov pri obrabotke davleniem [Resistance to deformation and plasticity of metals treated by pressure]. Moscow, Maschinostroenie Publ., 2011, 244 p.
[10] Shlykov Yu.P., Ganin E.A., Tsarevskiy S.N. Kontaktnoe termicheskoe soprotivlenie [Thermal contact resistance]. Moscow, Energiya Publ., 1977. 328 p.
[11] Yang C., Persson B.N.J. Contact mechanics: contact area and interfacial separation from small contact to full contact. J. of Physics: Condensed Matter, 2008, vol. 20, no. 21, 215214, 13 p. DOI:10.1088/0953-8984/20/21/215214
Статья поступила в редакцию 08.07.2015
Мурашов Михаил Владимирович — канд. техн. наук, доцент кафедры "Защита информации" МГТУ им. Н.Э. Баумана.
МГТУ им. Н.Э. Баумана, Российская Федерация, 105005, Москва, 2-я Бауманская ул., д. 5.
Murashov M.V. — Cand. Sci. (Eng.), Assoc. Professor of Information Protection department, Bauman Moscow State Technical University.
Bauman Moscow State Technical University, 2-ya Baumanskaya ul. 5, Moscow, 105005 Russian Federation.
Просьба ссылаться на эту статью следующим образом:
Мурашов М.В. Выбор параметров вычислительных алгоритмов при решении задачи контактного деформирования шероховатых тел в ANSYS // Вестник МГТУ им. Н.Э. Баумана. Сер. Естественные науки. 2016. № 1. C. 111-121. DOI: 10.18698/1812-3368-2016-1-111-121
Please cite this article in English as:
Murashov M.V. Selection of parameters for computational algorithms while solving the problem of rough bodies contact with ANSYS. Vestn. Mosk. Gos. Tekh. Univ. im. N.E. Baumana, Estestv. Nauki [Herald of the Bauman Moscow State Tech. Univ., Nat. Sci.], 2016, no. 1,pp. 111-121. DOI: 10.18698/1812-3368-2016-1-111-121
