CC BY
45
Моделирование процессов трения на основе совмещенного дискретно-континуального подхода
С.Г. Псахье, А.Ю. Смолин, Ю.П. Стефанов, Иг.С. Коноваленко
Институт физики прочности и материаловедения СО РАН, Томск, 634021, Россия
На основе дискретно-континуального описания проведено численное моделирование поведения пар трения, состоящих как из одинаковых, так и из существенно разнородных материалов. Отмечено образование «квазижидкого» слоя в парах из разнородных материалов в две стадии, первая из которых характеризуется образованием «сгустков», перемежающихся разреженными областями. Периодически меняющееся давление изменяет эффективный коэффициент трения, но амплитуда колебаний слабо сказывается на его значении. Сравнение результатов моделирования, полученных дискретно-континуальным и дискретным методами, показывает, что комбинированный метод позволяет объединить преимущества континуального и дискретного описаний и избежать их основных недостатков.
Simulation of friction processes on the basis of a hybrid discrete-continuum approach
S.G. Psakhie, A.Yu. Smolin, Yu.P. Stefanov, and Ig.S. Konovalenko
Institute of Strength Physics and Materials Science SB RAS, Tomsk, 634021, Russia
Using the discrete-continuum description we numerically simulate the behavior of friction pairs consisting of similar and of essentially different materials. It is noticed that in pairs of different materials a “quasi-liquid” layer is formed in two stages. The first stage is characterized by the formation of dense and rarefied zones. Periodically varying pressure changes the effective friction coefficient, but the oscillation amplitude exerts little effect on its value. Comparison between the results of discrete-continuum and discrete simulation shows that the hybrid approach allows one to combine advantages of the continuum and discrete descriptions and to avoid their main disadvantages.
1. Введение
Как известно, в процессах трения и износа в достаточно узкой зоне контакта двух тел происходят интенсивные процессы деформации, генерации и накопления повреждений, а также перемешивания. Моделирование именно этих процессов вызывает наибольшие трудности в методах механики сплошных сред. В то же время, имеются методы моделирования, основанные на дискретном описании, в которых упомянутые процессы учитываются естественным образом. Наиболее подходящим для данного класса задач представляется активно развиваемый в последние годы метод подвижных клеточных автоматов (МСА) [1-6]. Данный метод построен на дискретном описании, но в своей основе содержит несколько иные уравнения движения, отличные от классических для метода частиц. Основным достоинством данного метода является возможность явным образом моделировать как формирование несплошностей различного типа (от генерации отдельных повреждений до распространения магистральных трещин), так и эффект перемешивания масс.
Поскольку вне зоны контакта материалы трущихся тел деформируются незначительно, то для целей данной работы был разработан подход, в котором моделируе-
мый объект состоит из областей двух типов, одна часть которых описывается как сплошная, а другая — как дискретная среды (рис. 1). При этом для моделирования деформации континуальной части использовался конечно-разностный метод решения динамических задач упругопластического деформирования сплошных сред, а для моделирования дискретной части — метод подвижных клеточных автоматов. Выбранные методы обладают особенностями, позволяющими совместить их без особых затруднений. Оба основаны на лагранжевом подходе, то есть рассматривают движение каждой материальной точки в абсолютном пространстве, в обоих
Рис. 1. Дискретно-континуальная модель пары трения
© Псахье С.Г., Смолин А.Ю., Стефанов Ю.П., Коноваленко Иг.С., 2005
методах конечные уравнения движения могут быть записаны для точечных масс через действующие на них напряжения или силы. Кроме того, в [4] было показано, что при стремлении характерного размера автомата к нулю формализм метода МСА позволяет перейти к классическим соотношениям механики сплошной среды. Результатом явилась методика, позволяющая объединить преимущества обоих методов для моделирования процессов трения, включая описание генерации и накопления повреждений в условиях локальных интенсивных деформаций, реализующихся при динамических воздействиях в зонах контакта трущихся поверхностей. При этом дискретное описание применяется к узкой зоне контакта двух тел, а остальная часть образцов рассматривается в рамках механики сплошной среды.
2. Совмещение дискретного и континуального методов моделирования
Подробно методика совмещения дискретного и континуального методов, а также оба метода в отдельности изложены в работах [7, 8]. Суть совмещения состоит в том, что между дискретной и континуальной областями определяется некоторая граница сопряжения и полагается, что она принадлежит обеим областям. При этом каждому узлу расчетной сетки, лежащему на границе, ставится в соответствие некоторый автомат (элемент метода МСА). В общем случае, между двумя соседними узлами сетки, лежащими на границе, могут располагаться несколько клеточных автоматов. Тогда размер автоматов должен быть кратным шагу сетки: А = Бп, где А — шаг сетки; D — размер автомата; п — целое число. Перемещение граничных узлов вместе с находящимися в них сопряженными автоматами осуществляется в сеточном методе. Координаты автоматов, располагающихся между двумя соседними граничными узлами, рассчитываются с использованием интерполяции перемещений соответствующих узлов сетки.
Обеспечение непрерывности движения на поверхности раздела производится на этапе расчета скорости граничных узлов в континуальной области. Независимо от трактовки и используемого похода — конечно-разностного или конечно-элементного — при построении разностного аналога уравнений их форма и смысл слагаемых остаются одинаковыми. В левой части стоит произведение массы на компоненту ускорения узла рас-
четной сетки, а справа — компонента равнодействующей силы, приложенной к узлу:
X Г, т у = X Fiy, (1)
т х:
I=1, N
1=1, N
К =-аХх АУ +°Ху Ах',
Г = а Уу Ахг -а Ху А/.
Для всех внутренних узлов континуальной области выражение для суммарной силы, действующей на центр шаблона, показанного на рис. 2, N = 4, запишется в виде:
XГх = -(ахху24 + ахху31 + ахху42 + ахху13 ) +
I=ЦУ
. ^1 л,! ^ II . ТП III . ТУ 1Уч
+ (ахух24 + ахух31 + ахух42 + аху х13 )5
X Г = (аI х 1 + ап х 11 + ашх 111 + аIУх 1У ) -
1аГу ~ уух24 уух31 уух42 °уух13 )
(2)
I=I,IУ
Рис. 2. Шаблон построения разностного аналога уравнения движения. Показана нумерация узлов и ячеек
(гЛ I , —II II , _Ш III , _ДУ 1У \
- (ахуу24 +ахуу31 +ахуу42 +ахуу13 )■
Для всех граничных узлов слагаемые от недостающих ячеек заменяются силами, действующими на соответствующие автоматы со стороны МСА-области. Расчет движения таких узлов осуществляется исходя из уравнения (1), куда в зависимости от геометрии поверхности раздела будут входить силы, полученные в МСА-области. После перемещения автомата, неразрывно связанного с узлом сетки, соответствующее воздействие из области сетки передается во внутреннюю область МСА.
Таким образом, схематически алгоритм совместного расчета выглядит следующим образом:
1. В континуальной области решается система уравнений движения, осуществляется расчет напряженно-деформированного состояния, рассчитываются скорости и координаты.
2. Перед расчетом скоростей узлов, лежащих на границе раздела, осуществляется вызов подпрограммы, реализующей метод МСА. В нее передаются координаты и скорости совмещенных узлов автоматов с предыдущего шага, а также шаг интегрирования.
3. Пользуясь полученными данными, осуществляется шаг (а в случае мелких автоматов несколько шагов) интегрирования метода МСА. В результате рассчитываются новые положения и скорости всех автоматов, в том числе и силы, действующие на автоматы, совмещенные с узлами сетки.
4. Новые данные о граничных автоматах (узлах) возвращаются в сеточный метод. После этого в нем осуществляется расчет согласованного движения граничных узлов-автоматов, определяются их новые координаты.
5. Осуществляется согласование величины нового шага интегрирования по времени.
Разработанная методика была протестирована на задачах о распространении упругих волн различных типов в обеих средах и их прохождении через границу совме-
щения. Было показано, что в случае одинаковых свойств дискретной и континуальной частей, на границе совмещения не происходит искажения и отражения упругих волн [7].
3. Результаты моделирования
К настоящему времени в Институте физики прочности и материаловедения СО РАН выполнен комплекс работ по моделированию методом МСА процессов трения и износа стальных образцов [9, 10]. При этом область исследования ограничивалась узкой зоной контакта. Удалось показать, что в этой зоне может образовываться «квазижидкий» слой, где происходит интенсивное перемешивание вещества, определены зависимости толщины этого слоя и коэффициента трения от параметров нагружения.
В данной работе аналогичные расчеты были проведены с помощью совмещенной методики. Следует отметить, что требования к оперативной памяти существенно уменьшились и возросла скорость расчета. Это объясняется тем, что при расчете методом МСА модель включала 23470 автоматов, а в дискретно-континуальной модели использовалось 2860 автоматов и 1196 узлов сетки (рис. 3), а при одинаковом числе автоматов и узлов программа, реализующая сеточный метод, работает приблизительно в 5 раз быстрее. Значения толщины «квазижидкого» слоя и коэффициентов трения при расчете с использованием новой методики отличались от значений, полученных методом МСА, не более чем на 5 %. Эта разница объясняется различной реализацией граничных условий на верхней и нижней внешних гранях образцов.
Существенной особенностью перечисленных расчетов является то, что материалы обоих тел были одинаковыми (в данном случае сталь). Однако определенный интерес представляет случай, когда один из материа-
лов пары трения является существенно более жестким и прочным. Например, такие условия реализуются в искусственном тазобедренном суставе при трении головки эндопротеза (обычно делается из стали, титанового сплава или керамики) о вкладыш вертлужной впадины (полиэтилен) [11]. При моделировании подобных пар трения удобно одно из тел рассматривать абсолютно жестким. В данной работе были рассмотрены особенности приработки вкладыша на начальной стадии трения о головку эндопротеза. Общий вид расчетной области и параметры нагружения показаны на рис. 4 (дискретная область здесь и дальше показана в виде сетки межавтоматных связей). В данном случае вся нижняя часть образца, моделирующая поверхность головки, двигалась с заданной скоростью Vx, то есть считалась недеформируемой. Нас интересовало влияние прижатия (постоянного и периодического) в паре трения на процесс приработки. Размер автомата составлял 2.5 нм.
На первом этапе основными качествами материала вкладыша считаются: его относительная мягкость (малый модуль упругости по сравнению с материалом контртела) и склонность к необратимому деформированию. Поэтому в данной работе он моделировался ли-нейно-упрочняющейся пластической средой (модуль Юнга 12.6 ГПа, коэффициент Пуассона 0.3, предел текучести 120 МПа). В МСА-модели предел прочности составлял 150 МПа, а при достижении в несвязанной паре автоматов давления 152 МПа эта пара переходила в состояние связанности (то есть допускалась «перелинков-ка» автоматов в процессе перемешивания, как и в работах [9, 10]).
Результаты расчетов показаны на рис. 5. Следует отметить, что образование «квазижидкого» слоя в этом случае осуществляется в две стадии. На первой стадии затекание происходит не во все неровности жесткой поверхности (рис. 5, а). Причем за счет «перелинковки» автоматы из менее заполненных впадин могут переноситься в более заполненные, образуются «сгустки» из-
Рис. 3. Сравнение результатов моделирования методами МСА (а) и дискретно-континуальным (б)
Рис. 4. Общий вид дискретно-континуальной модели пары трения из разнородных материалов
Рис. 5. Ccto^ межавтоматных связей в моменты времени 0.108 (a) и 0.216 мкс (б)
ношенного материала, перемежающиеся разреженными областями. Однако по мере дальнейшего движения частицы мягкого материала заполняют все впадины (рис. 5, б) и распределяются достаточно равномерно вдоль контактной поверхности. Поверхность, разделяющая неповрежденный материал пластика от его изношенной части («квазижидкого» слоя), в данном случае является неровной, в отличие от трения тел из одинаковых материалов (рис. 3).
Анализ коэффициентов трения, полученных по результатам моделирования, показал, что они не зависят от амплитуды колебаний давления, прикладываемого к верхнему телу, но отличаются от значений, полученных при постоянном давлении (0.27 и 0.39 соответственно). Более интересным представляется вопрос о влиянии частоты колебаний на коэффициент трения, именно этой проблеме будут посвящены последующие исследования.
Таким образом, в данной работе на основе дискретно-континуального описания проведено численное моделирование поведения пар трения, состоящих как из одинаковых, так и из существенно разнородных материалов. Отмечено образование «квазижидкого» слоя в парах из разнородных материалов в две стадии, первая из которых характеризуется образованием «сгустков». При приложении периодически меняющегося давления амплитуда колебаний слабо сказывается на значении эффективного коэффициента трения. Показано, что совместное использование позволяет объединить преимущества континуального и дискретного описаний и избежать их основных недостатков.
Работа выполнена при финансовой поддержке Программы фундаментальных исследований Отделения энергетики, машиностроения, механики и процессов управления РАН № 3.12 (проект № 5) и РФФИ (проект № 05-08-33530-а).
Литература
1. Пcaxъe С.Г., Хори Я., Ko_ростелев С.Ю., Смолин И.Ю., Дмитриев A.И., Шилъко E.B., Aлeкceeв СМ. Mетод подвижных клеточных автоматов как инструмент для моделирования в рамках физической мезомеханики // Изв. вузов. Физика. - 1995. - Т. 38. -М 11.- C. 58-69.
2. Пcaxъe С.Г., Ko_ростелев С.Ю., Смолин A.Ю., Дмитриев AM., Шилъко E.B., MouceeнкoД.Д., Tamapuнцeв E.M., Aлeкceeв С.B. Mc-тод подвижных клеточных автоматов как инструмент физической мезомеханики материалов // Физ. мезомех. - 1998. - Т. 1. - М 1. -C. 95-108.
3. Пcaxъe С.Г., Дмитриев A.И., Шилъко E.B., Смолин A.Ю., Ko_росте-лев С.Ю. Mетод подвижных клеточных автоматов как новое направление дискретной вычислительной механики. I. Теоретическое описание // Физ. мезомех. - 2000. - Т. 3. - М 2. - C. 5-15.
4. Пcaxъe С.Г., Чертов M.A., Шилъко E.B. Интерпретация параметров
метода подвижных клеточных автоматов на основе перехода к континуальному описанию // Физ. мезомех. - 2000. - Т. 3. - М 3. -C. 93-96.
5. Popov V.L., Psakhie S.G., Dmitriev A.I., Shilko E. V. Quasi-fluid nanolayers at the interface between rubbing bodies: simulation by movable cellular automata // Wear. - 2003. - V. 254. - P. 901-906.
6. Попов B.Л., Пcaxъe С.Г., Жерве A., Kepвaлъд Б., Шилъко E.B., Дмитриев A-И. Износ в двигателях внутреннего сгорания: эксперимент и моделирование методом подвижных клеточных автоматов // Физ. мезомех. - 2001. - Т. 4. - М 4. - C. 73-83.
7. Пcaxъe С.Г., Смолин A.Ю., Cmeфaнoв Ю.П., Maкapoв KB., Шилъко E.B., Чертов M.A., Eвmушeнкo EM. Mоделирование поведения сложных сред на основе комбинированного дискретно-континуального подхода // Физ. мезомех. - 2003. - Т. 6. - М 6. - C. 11-21.
8. Смолин A.Ю., Cmeфaнoв Ю.П., Пcaxъe С.Г. ^вместное использование дискретного и континуального методов для моделирования деформации и разрушения в области контактного взаимодействия // Физ. мезомех. - 2004. - Т. 7. - ^ец. выпуск. - Ч. 1. - C. 70-73.
9. Попов B.Л., Пcaxъe С.Г., Шилъко E.B., Дмитриев A.И., K^me K., Бухер Ф., Эртц M. Исследование зависимости коэффициента трения в системе «рельс - колесо» как функции параметров материала и нагружения // Физ. мезомех. - 2002. - Т. 5. - М 3. - C. 17-25.
10. Popov V.L., Psakhie S.G., Dmitriev A.I., Shilko E. V. Mikro- und Na-noskala-Simulationen von Reibung und Verschleifi in Verbrennungs-motoren und in Rad-Schiene-Kontakten // Tribologie-Fachtagung 2002, Tribologische Systeme, Sept., Gottingen, Deutschland. - 2002. -33/1-33/7.
11. Бегун П.И., Шукейло Ю.A. Биомеханика. - CM.: Политехника, 2000. - 463 с.
