CC BY
88
УДК 53.072, 531.43, 539.62
Многоуровневое моделирование процессов трения и износа на основе численных методов дискретной механики и феноменологической теории
А.И. Дмитриев, А.Ю. Смолин, В.Л. Попов1, С.Г. Псахье
Институт физики прочности и материаловедения СО РАН, Томск, 634021, Россия Томский государственный университет, Томск, 634050, Россия 1 Берлинский технический университет, Берлин, D-10623, Германия
Предложена методика многоуровневого моделирования процессов трения и износа, объединяющая в рамках единой концепции численные расчеты методом подвижных клеточных автоматов на масштабах локальных контактов и макроскопическую феноменологическую теорию трения. При этом для моделирования на мезоскопическом уровне использована методика совмещения метода подвижных клеточных автоматов с численными методами механики сплошных сред. Приведены результаты моделирования локальных контактов в парах сопряжения гетерогенных фрикционных материалов. Проанализировано влияние исходной микрошероховатости на поведение квазижидкого слоя трения в пластичных материалах.
В рамках макроскопической феноменологической модели трения процессы выкрашивания, повторного приваривания частиц износа, а также процессы массопереноса от одного из контактирующих тел к другому описываются в форме стохастических дифференциальных уравнений. Параметры феноменологических уравнений могут быть найдены только в рамках микроскопических моделей процессов, происходящих в трибологическом контакте. В данной работе для этого предлагается использовать результаты моделирования методом подвижных клеточных автоматов. Согласно разработанной методике было найдено значение коэффициента диффузии слоя трения. Результаты проведенного исследования позволяют идентифицировать параметры макроскопической модели на основе моделирования процессов на наноуровне и тем самым «перебросить мост» между нано- и макроуровнями.
Ключевые слова: трение и износ, моделирование, методы частиц, дискретно-континуальное описание, феноменология
Multilevel simulation of friction and wear using numerical methods of discrete mechanics and phenomenological theory
A.I. Dmitriev, A.Yu. Smolin, V.L. Popov1, and S.G. Psakhie
Institute of Strength Physics and Materials Science SB RAS, Tomsk, 634021, Russia Tomsk State University, Tomsk, 634050, Russia 1 Technische Universitat Berlin, Berlin, D-10623, Germany
The paper proposes a multilevel simulation procedure for friction and wear, in which numerical calculation by the movable cellular automata method on scales of local contacts is combined with the macroscopic phenomenological theory of friction within the common conception. Mesoscopic simulation is performed by a procedure combining the movable cellular automata method with numerical methods of continuum mechanics. The simulation results for local contacts in conjugate pairs of heterogeneous friction materials are given. The influence of the initial microroughness on the behavior of a quasi-liquid friction layer in plastic materials is analyzed.
The macroscopic phenomenological model of friction employs stochastic differential equations to describe spalling, repeated welding of wear particles, and mass-transfer from one contacting body to another. Phenomenological equation parameters can be defined only by microscopic models of processes occurring in the tribological contact. To do this, the present paper proposes to use simulation results obtained by the movable cellular automata method. The developed procedure have allowed us to find the diffusion coefficient for the friction layer. The investigation findings make it possible to identify macroscopic model parameters by simulation of nanoscopic processes and thus “to throw a bridge” between the nano- and macrolevels.
Keywords: friction and wear, simulation, methods of particles, discrete-continued description, phenomenology
1. Введение
Согласно основным положениям физической мезо-механики поверхностный слой в деформированном твердом теле является самостоятельным мезоскопичес-
© Дмитриев А.И., Смолин А.Ю., Попов В.Л., Псахье С.Г., 2008
ким структурным уровнем деформации, который во многом определяет механическое поведение нагруженного материала в целом [1, 2]. Это объясняется не только исходной дефектностью и анизотропией поверхности,
но и возможностью реализации в приповерхностной области различных, в том числе недислокационных, механизмов деформации материала [2-5].
Исследование особенностей механизмов деформирования и разрушения поверхностных слоев приобретает особое значение при изучении процессов, связанных с взаимодействием контактирующих поверхностей и, в частности, триботехнических задач [4, 6, 7]. Несмотря на постоянное развитие экспериментальных методов материаловедения [3, 7] зона трения и, особенно, зона фактического контакта, остаются труднодоступными для исследования непосредственно в процессе испытания. Как правило, основные данные могут быть получены только после завершения эксперимента или в какой-то промежуточный момент после остановки испытания. Сложность и многофакторность процессов, происходящих в зоне контакта, не позволяет получить однозначной связи между структурными, физико-механическими и триботехническими характеристиками материала. Для понимания закономерностей деформирования и разрушения поверхностного слоя при трении и износе эффективным представляется использование методов численного моделирования [8-11]. Результаты, полученные в процессе моделирования, могут составить основу для прогнозирования поведения материалов при контактных взаимодействиях и обеспечить в дальнейшем совершенствование свойств триботехнических материалов. В свою очередь, сложность проблемы теоретического описания состоит не только в разнообразии микромеханизмов трения и износа, но также и в многомасштабности задачи. Так, согласно исследованиям микроструктуры и химического состава приповерхностных слоев элементарные процессы износа происходят на масштабном уровне порядка 10-100 нм [10], в то время как макроскопические изменения топографии поверхности и транспорт частиц износа из области трибологического контакта происходят на масштабах порядка 1-10 мм. Оба эти масштаба являются существенными и определяющими для процессов износа. Тот факт, что они различаются на 5-6 порядков, делает в настоящее время невозможным описание процессов трения и износа в рамках единой модели, включающей как фундаментальные физические процессы на нано-метровом уровне, так и износ как таковой.
В рамках настоящей работы нами предлагается использовать многоуровневый подход для моделирования процессов деформации и формирования повреждений в поверхностном слое трибосопряженной пары на разных масштабах. Суть подхода состоит в использовании разномасштабных методик исследования в рамках единой структуры данных. Ниже приведены основные принципы развиваемых подходов и результаты численного моделирования контактных задач, полученные в рамках предлагаемых методик.
2. Особенности компьютерного моделирования задач трения и износа
Для численного моделирования процессов трения и износа на уровне пятен локальных контактов необходимо принять во внимание характерный масштаб, на котором в поверхностных слоях контактирующих тел развиваются механизмы, определяющие эти процессы. При этом надо учесть, что трение и износ сопровождаются нарушениями сплошности и интенсивным мас-сопереносом. По этой причине для изучения трибологических процессов нанометрового диапазона нами был использован численный метод дискретного подхода -метод подвижных клеточных автоматов (МСА). Формализм метода подробно описан в [12]. Основным достоинством данного метода является возможность явным образом моделировать как формирование несплошнос-тей различного типа (от генерации отдельных повреждений до распространения магистральных трещин), так и эффекты перемешивания масс [9, 10, 12, 13].
Очевидно, что выбранный пространственный масштаб исследований определяет величину шага интегрирования численной схемы, который находится в диапазоне 1-10 пс. Понятно, что такие малые пространственные и временные диапазоны не позволяют рассматривать методом подвижных клеточных автоматов протяженные объекты на длительных временных интервалах. Применение идеи разномасштабных методик исследования в рамках единой системы данных позволяет эффективно решить эту проблему.
В работе предлагаются две методики: 1) совмещение дискретного и континуального подходов описания моделируемой среды с разделением областей по интенсивности протекающих в них процессов деформации и разрушения и 2) совместное использование микроскопического моделирования с макроскопической феноменологической теорией трения, формулируемой в форме стохастических дифференциальных уравнений. Рассмотрим обе методики более подробно.
3. Дискретно-континуальное описание
Как было показано выше, непосредственная зона контакта при трении может быть успешно смоделирована с помощью метода подвижных клеточных автоматов. Однако при перечисленных преимуществах дискретного подхода у него есть и недостатки. В частности, использование малого размера автоматов для детального описания зоны контакта приводит к существенным требованиям к оперативной памяти компьютера. Способность методов частиц описывать интенсивное перемешивание вещества требует пересчета на каждом шаге всего окружения для каждой частицы и, значит, большого времени счета. Поэтому в качестве первого подхода многоуровневого моделирования предлагается исполь-
зовать процедуру совмещения континуального и дискретного описания моделируемой среды [14]. Континуальное представление применяется для описания тех участков материала, которые деформируются преимущественно в упругой области и вероятность формирования различного рода повреждений в них очень мала, а метод подвижных клеточных автоматов — для описания поведения областей, подверженных наибольшим воздействиям с возможностью развития процессов разрушения и перемешивания материала. При этом в области континуального описания можно использовать расчетные ячейки в несколько раз большие, чем размер подвижных клеточных автоматов. Использование метода конечных элементов позволяет использовать адаптивную сетку, в которой элементы, находящиеся далеко от границы совмещения, могут быть еще больше.
Рассмотрим некоторые результаты, полученные в рамках предложенного подхода. В качестве континуального метода использовался конечно-разностный метод решения динамических задач механики деформируемого твердого тела. Сам метод и подробная процедура его совмещения с методом подвижных клеточных автоматов представлены в [14].
На примере пары трения «тормозной диск - колодка» были смоделированы следующие экспериментально наблюдаемые контактные ситуации: взаимодействие двух хрупких и взаимодействие двух пластичных модельных материалов. Схема контакта, а также области, моделируемые различными методами, показаны на рис. 1, а. Первоначально автоматы, соответствующие разным телам, располагались на некотором расстоянии. Затем ко всем автоматам нижнего блока прикладывалась постоянная горизонтальная скорость V, линейно увеличивающаяся от 0 до 10 м/с (рис. 1, б). После достижения максимального значения скорости на автоматы верхнего слоя верхнего блока начинала действовать постоянная сила, имитируя сжатие. Максимальное значение силы соответствовало локальному давлению Р в диапазоне 20-100МПа. В горизонтальном направлении задавались периодические граничные условия. Слои металлической подложки (верхний для колодки и нижний для диска), подверженные преимущественно деформациям в упругой области, моделировались с помощью метода континуальной механики. При этом для строгого выполнения условий непрерывности движения на границах совмещения верхняя часть верхнего МСА-блока и нижняя часть нижнего МСА-блока задавались автоматами с параметрами соответствующей стали: ферритная сталь для диска и перлитная — для колодки.
Следует отметить, что в случае контакта пластически деформируемых материалов наряду с возможностью разрыва межавтоматных связей в модели задавалась возможность формирования новых связей между отдельными автоматами, принадлежащими как одному, так и разным телам. Этот переход может быть интерпрети-
Сетка
Рис. 1. Начальная структура (а) и схема нагружения моделируемой пары трения (б)
рован как развитие процессов наносварки в пятнах контакта. В качестве критерия переключения пары «несвязанных» автоматов в состояние «связанные» использовалось достижение локальным напряжением в паре величины предела упругости модельного материала.
Результаты моделирования показывают, что тип контакта определяет структуру граничного слоя. Так, для контакта модельных хрупких материалов это преимущественно независимые частицы износа, механически перемешиваемые во время относительного движения (рис. 2, а). Топология поверхностей, формируемых в процессе относительного движения, состоит из множества оборванных цепочек автоматов и имеет выраженный фрактальный характер. При моделировании
контакта двух сталей граничный слой остается связанным с обоими блоками, однако его структура характеризуется большим количеством дефектов (рис. 2, б). Такое отличие объясняется не только свойствами материалов, но и учетом возможности восстановления связей для контакта сталей.
Различие наблюдается не только в структуре формируемого граничного слоя, но и в величине силы сопротивления относительному проскальзыванию, т.е. в рассчитываемом значении коэффициента трения для рассматриваемых случаев. Характерно, что для контакта хрупких материалов коэффициент трения находится в диапазоне 0.3—0.5 в зависимости от наличия в них дополнительных включений, а для контакта пластически деформируемых материалов его значение достигает 0.8, что находится в хорошем соответствии с экспериментальными данными.
Далее рассмотрено влияние параметров исходного профиля контактирующих поверхностей и скорости относительного движения на структуру граничного слоя и коэффициент трения. Кроме того, рассматривался случай наличия в обоих телах неоднородной периодической структуры из прочных гексагональных включений и относительно податливой прослойки. На рис. 3 представлены МСА-области некоторых моделируемых образцов в различные моменты времени.
Расчеты показывают, что во всех рассмотренных случаях коэффициенты трения практически совпадают, а толщина граничного слоя отличается не более чем на 50 %. В период формирования этого слоя поверхности образцов образуют достаточно сложный рельеф за счет процессов выкрашивания частиц износа и их приваривания (наносварки). При этом в однородных образцах прослеживается тенденция к образованию макрорельефа с характерным размером, равным ширине расчетной области (вдоль направления, в котором заданы периодические граничные условия). Этот факт подтверждается наличием соответствующих частот в спектре упругих колебаний, генерируемых в процессе трения. В неоднородных материалах в спектре регистраций интенсивности напряжений присутствуют частоты, соответствующие смещению взаимодействующих тел на характерный размер неоднородностей.
Таким образом, результаты компьютерных исследований с помощью метода подвижных клеточных автоматов (см. также [9, 10, 13]) показывают формирование в области локального контакта так называемого квази-жидкого слоя, в котором протекают интенсивные динамические процессы, ответственные за перенос масс контактирующих материалов и изменение профиля контактирующих поверхностей. Очевидно, что аналогичные процессы могут быть рассмотрены и для больших
Рис. 3. Конфигурации межавтоматных связей для пар трения различной структуры при скорости относительного движения 2 м/с в моменты времени t = 0 (а), 1 (б) и 4 мкс (в)
систем, для которых непосредственное моделирование изменения профиля поверхности и расчета поверхностного спектра мощности с помощью метода подвижных клеточных автоматов становится очень трудоемкой задачей. С этой целью при переходе на большие масштабы предлагается использование феноменологической модели трения [15]. В рамках данной модели предлагается условно разделить динамические процессы, протекающие в квазижидком слое, на те, что ответственны за стохастический массоперенос материала вдоль поверхности контакта, и те, которые приводят к изменению профиля поверхности.
4. Методика совмещения с феноменологической моделью трения
4.1. Формализм макроскопической феноменологической модели трения
В первом приближении механизмы, ответственные в процессе трения и износа за транспорт материала сопряженных поверхностей, могут быть описаны случайными «прыжками» с изменением функции высоты ^х, г) конкретных рассматриваемых дискретных точек поверхности х.. В этом случае функция высоты для точки х. выбранной поверхности трения в момент времени г может быть записана как к(х{, г) = к, (г). Ее значение со временем меняется стохастически, что может быть интерпретировано как результат пластического деформирования. Таким образом, если высоты в точках (;' + 1) или (;' - 1) возрастают на величину Дh, то согласно закону сохранения массы это приводит к тому, что высота в точке i уменьшается на то же значение, например как результат соударения с частицей износа. Правило перехода для описанного процесса в вероятностной форме может быть сформулировано в следующем виде:
.>...........к. - к —
И 1 = и0— + ------—, (1)
01 Д х
где j = i + 1 или j = i - 1. Последний множитель в (1) содержит зависимость вероятности перемещения массы от величины градиента высот. Этот дискретный процесс может быть представлен в рамках основного уравнения, описывающего вероятность распределения функции высоты. С использованием формализма Фоккера-План-ка процесс переноса масс может быть записан в форме стохастического дифференциального уравнения следующего вида [15, 16]:
д д ^ д
—к(x, 0 = Г>1 — к(x, 0 + — —(х г) (2)
дг дх дх
Д — коэффициент «диффузии», который связан со скоростью массопереноса как
и1Дх2 Дх^° > D1, (3)
а стохастический поток —8 имеет корреляционную
функцию, соответствующую так называемому «белому шуму», т.е. равновероятностному распределению:
— (х, 0—* (хг, 0) = 2DoS(x - х) 8(г - г'), (4)
D0 определяется как
и0Дх4 Дх>х > D0. (5)
В (4) треугольные скобки означают усреднение по всему ансамблю.
Следует отметить, что описываемая таким образом «диффузия» учитывает только механически активируемый перенос масс и не учитывает термические воздействия, что, безусловно, для целого ряда задач может иметь первостепенное значение.
Следующим процессом, рассматриваемым в рамках предлагаемой макроскопической феноменологической модели трения, является учет изменения профиля сопряженных поверхностей, что может быть представлено в виде случайных событий: отрыва частиц, принадлежащих поверхности (микроизнос), с выносом их в слой трения и случайным присоединением частицы износа к выделенной поверхности контактирующего материала (микросварка, адгезия и др.). Эти процессы описываются стохастической функцией g8((х, г), симметричной относительно событий отрыва и присоединения частиц. Это предположение эквивалентно тому, что квазижид-кий слой имеет постоянную толщину и среднее от gst равно нулю. Кроме того, делается предположение, что стохастические воздействия со стороны слоя трения являются нерегулярными, нескоррелированными на малой длине рассматриваемого временного интервала, что для протяженной системы означает, что корреляционная функция соответствует распределению вида «белого шума»:
^(х г^ = 0 (6)
(gst(x, г)gst(x,, г')) = 2^(х-х)8(г-г").
Это выражение также может быть выведено с использованием формализма Фоккера-Планка. Корреляционный коэффициент G является мерой силы ударов, приводящих к процессам наноизноса и наносварки на различных участках поверхности, и может быть интерпретирован как скорость «прыжков»:
Дх2 Дх>0 > G. (7)
Следует отметить, что коэффициент G связан с коэффициентом «диффузии» D, который в свою очередь характеризует вертикальные флуктуации положения ква-зижидкого слоя. Поскольку коэффициент G связан с вероятностью прыжков в одной конкретной точке, а D характеризует флуктуации положения всего слоя трения с конкретной его толщиной, то связь между G и D может быть записана как
D = , (8)
л/я
где п — толщина квазижидкого слоя.
Учитывая вышесказанное, полное уравнение для предлагаемой феноменологической модели может быть
представлено в виде стохастического дифференциального уравнения, описывающего динамическое изменение топологии поверхности трения с феноменологическими константами D0, D1 и О, определяемыми согласно выражениям (5), (3) и (7) соответственно: д
— к (х, г) = £й(X, г) +
дг
д 2 д
+ А ТГк(х>г) + ^ —81(х>г). (9)
дх 2 дх
Это уравнение позволяет проследить за динамикой развития поверхности трения. Как показывает эксперимент, трение и износ приводят к формированию фрактальных поверхностей [17]. Параметры феноменологических уравнений могут быть найдены только в рамках микроскопической модели процессов, происходящих в трибологическом контакте. Случайное блуждание положения поверхности износа, которое может быть описано эффективным коэффициентом «диффузии» является одним из важнейших параметров названной феноменологической модели. Далее приведены результаты исследования этих случайных блужданий и методика оценки изменения профиля поверхности трибосопря-жения исходя из микроскопической модели трения, основанной на методе подвижных клеточных автоматов. Результаты этих исследований позволяют идентифицировать параметры макроскопической модели на основе моделирования процессов на наноуровне и тем самым «перебросить мост» между нано- и макроуровнями.
4.2. Оценка коэффициента «диффузии» слоя трения на основе моделирования методом подвижных клеточных автоматов
В расчетах использовалась рассмотренная ранее схема нагружения (рис. 1), имитирующая локальный контакт в паре трибосопряжения. Для имитации протяженности моделируемых блоков в вертикальном направлении использовалась процедура диссипации энергии по вязкоупругому закону для автоматов граничных слоев. В горизонтальном направлении моделировались периодические граничные условия, имитируя тем самым
замкнутую систему трения. Начальная шероховатость моделируемых поверхностей задавалась в явном виде, однако, как показано выше, ее изменение в пределах глубины слоя шероховатостей максимального размера не оказывает существенного влияния на конечные результаты моделирования. Выбор размера автомата был обусловлен экспериментально наблюдаемым характерным размером фрагментированной структуры материалов трения, используемых для производства автомобильных тормозных систем. Размеры таких частиц составляют десятки нанометров [10], поэтому размер автоматов моделируемого образца d был равен 10 нм. Рассматривался сухой контакт модельных материалов со свойствами перлитной стали. Параметры модельных материалов и вид функции отклика автомата приведены в табл. 1. Размер автомата и упругие свойства модельных материалов определили шаг интегрирования численной схемы Дг, который составил 2.5 • 10-13 с. Следует отметить, что при моделировании наряду с разрывом межавтоматных связей задавалась возможность восстановления сплошности материала между независимыми элементами среды, что имитировало процессы плавления, наносварки и пр. Критерием разрыва межавтомат-ной связи служило достижение рассчитываемой интенсивностью напряжения в паре значения, равного прочности модельного материала. В качестве критерия для восстановления разорванной связи вступало достижение напряжением на контакте двух автоматов величины (О8 -ау2)/2 (табл. 1).
Анализ структуры поверхности взаимодействующих тел показал, что в процессе их относительного движения формируется слой, в котором происходят интенсивные процессы разрушения и сглаживания шероховатостей, разрыва существующих и создания новых связей между автоматами, а также перемешивание частиц материалов контактирующих тел. Кроме того, 7-координата слоя трения не является строго фиксированной величиной и может меняться в процессе взаимодействия (рис. 4, б).
Для оценки эффективного коэффициента «диффузии» слоя трения в процессе моделирования вычислялось среднее значение 7-координаты всех автоматов
Параметры модельного материала и вид функции отклика
Таблица 1
Модуль Юнга Е, ГПа
Коэффициент Пуассона V
Предел упругости а у1, МПа
Предел текучести а у2, МПа
Прочность а 8, МПа
Деформация при пределе текучести в у2
Предельная деформация на разрыв в 8
206
0.28
520
800
920
0.04
0.106
Рис. 4. Фрагмент структуры межавтоматных связей на этапе установившегося движения (а) и изменение положения квазижидкого слоя в процессе трения (б)
интерфейсной области, не имеющих связи ни с одним из блоков. Среднеквадратичное смещение 7-координаты слоя от первоначального значения вычислялось по формуле [18]:
(л (г )2) = NЕу (г) - у (г0)]2, (10)
где за начало отсчета выбирался момент времени, соответствующий установлению в системе динамического равновесия; N — число расчетов в серии; У. (г) и У (г 0) — текущее и начальное значения 7-координаты слоя трения для /‘-го расчета. Результат усреднялся по серии из 32 расчетов. Расчеты в серии отличались случайным разбросом параметров модельного материала для автоматов моделируемого образца по нормальному закону в диапазоне 10 % при идентичных условиях нагружения.
Согласно полученным данным в ходе относительного движения контактирующих поверхностей наблюдается дрейф граничного слоя трения с показателем степени около 0.9 (рис. 5, а). Для сравнения на рис. 5, а также
приведена линейная зависимость. При этом наблюдаются участки «быстрой диффузии» с показателем степени до 2 и участки относительного «покоя». Анализ структуры межавтоматных связей показал, что участки быстрого дрейфа связаны с вовлечением новых неповрежденных слоев автоматов или их частей.
Коэффициент «диффузии» D может быть найден из соотношения:
lim ^ R (t) 2 ^ = const + 2 Dt. (11)
Для рассматриваемой системы длиной 600 нм рассчитанное значение коэффициента «диффузии» D = 2.1 х х10-11 м2/с (рис. 5, б).
Анализ динамики граничного слоя трения показал, что толщина слоя медленно растет со временем (рис. 6, а). Изменение толщины слоя при установившемся режиме трения хорошо описывается степенным законом с показателем степени 0.25. Аналогичный нелинейный характер кривой прослеживается и для зависимости изменения числа частиц износа от времени
Рис. 5. а — изменение среднеквадратичного смещения слоя трения со временем, среднеквадратичное смещение приведено в единицах d2, время указано в единицах интервала записи конфигурации клеточных автоматов п, равного 104 • Дг = 2.5 нс; б — изменение эффективного коэффициента диффузии со временем
f"1 record f"1 record
Рис. 6. Изменение толщины слоя трения h в единицах диаметра автомата d (а) и изменение числа частиц износа nwear (б) со временем в логарифмических единицах
(рис. 6, б). Число частиц износа рассчитывалось как сумма автоматов, определяющих положение слоя трения, т.е. несвязанных ни с верхним, ни с нижним блоками.
4.3. Анализ изменения топологии поверхности трения
Для анализа топологии сформировавшихся поверхностей построим кривую, описывающую все пары автоматов, находящиеся в состоянии «связанные» с одной из контактирующих поверхностей. На рис. 7 такая кривая изображена для нижней поверхности. С помощью преобразования Фурье соответствующую этой кривой функцию можно представить суммой гармонических составляющих (синусоид и косинусоид с различными амплитудами А, периодами Т или частотами ш). Поскольку в методе подвижных клеточных автоматов материал, а следовательно, и его поверхность представлены в виде набора дискретных элементов, воспользуемся дискретным преобразованием Фурье [19].
Схема дискретного представления искомой функции показана на рис. 8. Величина шага дискретизации 8 соответствует радиусу автомата (минимальное расстояние между вертикальными проекциями центров автоматов в плотной упаковке, расположенных друг над другом). Минимальное расстояние между горизонтальными проекциями центров автоматов соответствует шагу
шероховатости рассматриваемой поверхности. Считаем, что на протяжении шага дискретизации значение функции постоянно и равно значению координаты автомата 7, находящемуся в данном интервале. Согласно преобразованию Фурье выражение для спектра мощности записывается в виде:
с 2( Кп) = а 2( Кп) + ь 2( Кп), (12)
где Кп = 2пп//, а целые п е (0,..., М), М = 8//. Выражения для а2(Кп ) и Ь 2(Кп ) в случае, когда f постоянна на длине шага дискретизации, могут быть записаны в виде сумм:
2 1
- J f(x)sin( Knx)<h-l 0
M mSx
^ fm J sin(KnX)dx
(m-1) Sx
m=1
l
b (Kn) = - J f (x)cos( KnX )dx =
(13)
mSx
И
Е, ^ / С08( Кпх ^
т=1 (
(т-1) 8х
Поскольку интегралы в правой части выражения (13) могут быть вычислены явно, перепишем значения для функций а (Кп ) и Ь (Кп ) как
Рис. 7. Пример функции, описывающей профиль нижней поверхности пары трибосопряжения
Рис. 8. Алгоритм дискретизации функции, описывающей профиль поверхности, с учетом шероховатостей, соответствующих размеру автомата
-5.0
-5.5
-6.0
...-6.5
<5 -7.0
О)
-7.5
-8.0
-8.5
-9.0
\а_
■ 1 / 2
/ ’
!к
ш ч
ч
и
1.5 2.0 2.5 3.0 3.5 4.0
<|д Ю
<|д к)
Рис. 9. Спектральные плотности поверхностных профилей: для верхней (1) и нижней (2) поверхностей при внешних условиях нагружения Fy = 255 МПа, ¥х = 3 м/с (а); для разных давлений: 127 (1), 255 (2), 383 МПа (3) при ¥х = 5 м/с (б); для разных скоростей: 3 (1), 4 (2), 5 м/с (3) при Еу = 383 МПа (в)
а (Кп) =
Ь (Кп) =
И
V*
т=1
И
Е/т
т=1
Кп
С08
(Кпх )|
т8х
(т-1) 8х
^п (Кпх)|
Кп
т8х
(т-1)8х
(14)
с (Кп) =
1 т8х
-еХр(гКпх)1(т-1)&
(15)
Согласно (12) функция с (Кп ) запишется в виде:
оГ И (
2 Е
т=1 1Кп п
где обе фазы данного комплексного коэффициента меняются стохастически.
Основная информация о поверхности может быть получена из спектральной плотности, которая рассчитывается как
12 1
£ (Кп ) = |с (Кп )|:
2п
(16)
На рис. 9, а показаны спектральные плотности верхней и нижней поверхностей для случая V = 3 м/с и силы Fy, соответствующей давлению 255 МПа в логарифмических координатах. Во избежание флуктуационных отклонений приведенные значения для ^(£ (Кп )) усреднялись по пяти профилям поверхности в последовательные моменты времени. Спектральная плотность исследуемых поверхностей показывает типичную зависимость волнового вектора с углом наклона кривой, приблизительно равным двойке, что справедливо для многих поверхностей фрактальной шероховатости.
На рис. 9, б и в показаны спектральные плотности для поверхностей в случае варьирования значения приложенного давления и скорости относительного движения сопряженных поверхностей. Хорошо видно, что с ростом давления кривые располагаются выше. Это означает, что большее давление приводит к формиро-
ванию более шероховатой поверхности. В то же время, с ростом Ух эта разница уменьшается, что говорит о том, что с повышением относительной скорости проскальзывания влияние наношероховатостей поверхности повышается вне зависимости от приложенного давления, и профили поверхностей становятся подобными.
Очевидно, что на основе данных численного эксперимента, выполненного с помощью метода подвижных клеточных автоматов для протяженных сопряженных поверхностей, можно провести исследование динамики изменения профиля поверхности, характеризующейся скоростью смещения точки излома кривой спектра мощности. Это, в свою очередь, определяет интенсивность макроскопических процессов, ответственных за изменение топологии сопряженных поверхностей (один из механизмов феноменологической модели трения) и позволяет на основе микроскопических исследований оценить величину коэффициента О.
Таким образом, полученные результаты позволяют идентифицировать параметры феноменологической макроскопической модели трения, а также с новых позиций рассматривать динамику процессов, протекающих на поверхности контактирующих тел при трении и износе. Следует отметить, что использование периодических граничных условий накладывает определенные ограничения на оценку макроскопических параметров феноменологической модели. В общем случае толщина слоя трения зависит от вертикальной подвижности квази-жидкого слоя, что определяется коэффициентами D1 и О. В условиях использования замкнутой системы трения, т.е. периодических граничных условий, вертикальные перемещения слоя трения будут определяться только параметром О.
5. Заключение
Представленные результаты демонстрируют возможность многоуровневого описания трибологических задач в рамках единой структуры данных. Такой многомасштабный подход позволяет более детально и всесторонне исследовать динамику процессов, протекающих на поверхности пятен контакта при трении. В частности, при решении контактных задач это позволяет учитывать механизмы микроскопического масштаба, ответственные за формирование граничного слоя трения, в котором интенсивно протекают процессы массопереноса, деформации, пластического течения, разрушения и сглаживания шероховатостей, и исследовать последствия этих процессов на более длительных временных интервалах и для макроскопически протяженных объектов, что имеет существенно большее значение с практической точки зрения. Кроме того, использование многоуровневого подхода позволяет явно моделировать поведение всего объема контактирующих тел, что является принципиально важным для корректного решения динамических задач трения, поскольку дает возможность явно учитывать глубину проникновения периодических воздействий со стороны поверхности с длинами волн от размера одной микронеровности до ширины моделируемой области контакта. Именно возможность учета контактирующих материалов как целого позволяет учесть влияние всего спектра колебаний — от элементарной шероховатости до размера протяженных объектов пары трения.
Важным результатом работы является также вычисление двух параметров, характеризующих свойства ква-зижидкого слоя трения. Это толщина слоя частиц износа и эффективный коэффициент «диффузии», описывающий интенсивность стохастических процессов массопе-реноса в зоне трения. Оба эти параметра позволяют рассчитать скорость износа в макроскопической трибологической системе с учетом конкретных размеров и геометрии системы. Вычисление этих параметров является нетривиальной задачей, поскольку она должна решаться на мезоскопическом уровне (характерная длина 10 нм) и учитывать не только процессы пластической деформации и разрушения, но и массопереноса и перемешивания. Несомненно, что для более глубокого понимания и детального анализа данных вопросов необходимо проведение дополнительных теоретических и экспериментальных исследований на различных масштабных уровнях.
Работа выполнена при финансовой поддержке проекта № 4.13.1 Программы специализированных отделе-
ний РАН, гранта Немецкой службы академических обменов DAAD Ref. 325 и гранта РФФИ № 07-08-00192-а.
Литература
1. Панин В.Е. Физическая мезомеханика поверхностных слоев твердых тел // Физ. мезомех. - 1999. - Т. 2. - № 6. - С. 5-24.
2. Панин В.Е., Фомин В.М., Титов В.М. Физические принципы мезо-
механики поверхностных слоев и внутренних границ раздела в деформируемом твердом теле// Физ. мезомех. - 2003. - Т. 6. -№2. - С. 5-14.
3. Veprek S., Jilek M. Superhard nanocomposite coatings. From basic science toward industrialization // Pure Appl. Chem. - 2002. - V. 74. -No. 3. - P. 475-481.
4. Крагелъский И.В. Трение и износ. - М.: Машгиз, 1962. - 384 с.
5. Rozman M.G., Urbakh M., Klafter J. Stick-slip dynamics of interfacial
friction // Physica A. - 1998. - V. 249. - P. 184-189.
6. Браун Э.Д., Буше Н.А., Буяновский И.А. и др. Основы трибологии (трение, износ, смазка): Учебник для технических вузов / Под ред. А.В. Чичинадзе. - М.: Наука и техника, 1995. - 778 с.
7. Schofer J., Santner E. Quantitative wear analysis using atomic force microscopy // Wear. - 1998. - V 222. - P. 74-83.
8. ГорячеваИ.Г. Механика фрикционного взаимодействия. - М.: Нау-
ка, 2001. - 478 с.
9. Dmitriev A.I., Popov V.L., Psakhie S.G. Simulation of surface topography with the method of movable cellular automata // Tribology Int. - 2006. - V. 39. - No. 5. - P. 444-449.
10. Osterle W., Dmitriev A., Klofi H., Urban I. Towards a better understanding of brake friction materials // Wear. - 2007. - V. 263. - P. 11891201.
11. Kemmer H.A. Investigation of the Friction Behavior of Automotive Brakes through Experiments and Tribological Modeling / Ph.D. Thesis. - Paderborn: University of Paderborn, Bosch publ., 2002. - 280 p.
12. Псахъе С.Г, Остермайер Г.-П., Дмитриев А.И., Шилъко Е.В., Смолин А.Ю., Коростелев С.Ю. Метод подвижных клеточных автоматов как новое направление дискретной вычислительной механики. I. Теоретическое описание // Физ. мезомех. - 2000. -Т. 3. - №2. - С. 5-13.
13. Popov V.L., Psakhie S.G., Shilko E. V, Dmitriev A.I. Quasi-fluid nanolayers at the interface between rubbing bodies: simulation by movable cellular automata // Wear. - 2003. - V. 254. - No. 9. - P. 901-906.
14. Псахъе С.Г., Смолин А.Ю., Стефанов Ю.П., МакаровП.В., Шилъко Е.В., ЧертовМ.А., ЕвтушенкоЕ.П. Моделирование поведения сложных сред на основе комбинированного дискретно-континуального подхода // Физ. мезомех. - 2003. - Т. 6. - № 6. - С. 1121.
15. Schargott M., Popov V.L. Diffusion as a model of formation and development of surface topography // Tribology Int. - 2006 - V. 39. -No. 5. - P. 431-436.
16. Gardiner C.W Handbook of Stochastic Methods for Physics, Chemistry, and Natural Sciences / Ed. by H. Haken. - Berlin-Heidelberg-New York-Tokyo: Springer, 1983. - 442 p.
17. Barabasi A.L., Stanley H.E. Fractal Concepts in Surface Growth. -Cambridge: Cambridge University Press, 1995. - 366 p.
18. Рит М. Наноконструирование в науке и технике. Введение в мир нанорасчета. - Москва-Ижевск: НИЦ «Регулярная и стохастическая динамика», 2005. - 160 с.
19. Поттер Д. Вычислительные методы в физике. - М.: Мир, 1975. -218 с.
Поступила в редакцию 13.12.2007 г.
Сведения об авторах
Дмитриев Андрей Иванович, д.ф.-м.н., ведущий научный сотрудник ИФПМ СО РАН, [email protected] Смолин Алексей Юрьевич, к.ф.-м.н., старший научный сотрудник ИФПМ СО РАН, [email protected] Попов Валентин Леонидович, д.ф.-м.н., профессор Берлинского технического университета, [email protected] Псахье Сергей Григорьевич, д.ф.-м.н., профессор, директор ИФПМ СО РАН, [email protected]
