Modeling of sliding contact of a system of asperities and coated viscoelastic half-space

1. Введение. Отказы техники часто связаны с неудовлетворительной работой узлов трения, изготовленных из высокоэластичных материалов, в частности уплотнений. Причиной использования высокоэластичных материалов, чаще всего это резины, является необходимость герметизации и поглощения вибраций различной природы (в том числе возникающих при трении). Однако эластомерные материалы без специальных силовых элементов неработоспособны при высоких давлениях и скоростях скольжения. Нанесение покрытий является одним из основных методов модификации для обеспечения антифрикционных свойств и износостойкости, при этом покрытия должны быть относительно тонкими для сохранения герметизирующих свойств уплотнения. Толщина покрытия может быть сравнима с размерами неровностей поверхностей взаимодействующих тел, поэтому задача множественного контакта вязкоупругих тел с покрытием в условиях скольжения является практически важной и актуальной.

Для вязкоупругого материала деформационный механизм трения напрямую связан с реальной шероховатостью поверхности, которую трудно математически описать, поскольку она имеет сложную многомасштабную природу. В последнее время появились значимые работы для описания самоподобных шероховатых поверхностей [1––4], которые обычно образуются при случайном характере формирования их геометрии. Многие исследователи используют подход, предложенный Перссоном в работе [3], для математического описания контакта самоподобной шероховатой поверхности. Численные подходы для моделирования скольжения резины по шероховатой поверхности были разработаны в работах [5––9]. В случае, когда геометрия поверхности возникает в результате шлифования или точения модели, ее нельзя описать нормальным распределением Гаусса, поэтому модель контакта, опирающуюся на предположение об около-фрактальной природе поверхности сложно брать в расчет [10]. Для моделирования деформационной составляющей силы трения используются различные подходы. Скольжение стандартных линейных вязкоупругих тел по шероховатым поверхностям рассматривалось в работах [2, 3]. Полученные результаты хорошо согласуются с экспериментальными данными. В работе [11] рассматривалось скольжение шероховатого индентора по вязкоупругому основанию, описываемому моделью Кельвина. Поверхность состояла из самоаффинных прямоугольных выступов. Исследована деформационная составляющая силы трения в зависимости от нагрузки, скорости скольжения и геометрических характеристик шероховатости. Хотя фрактальный подход эффективен для моделирования рельефа реальных поверхностей, рассматриваемые в упомянутых работах модели вязкоупругих материалов достаточно упрощены.

Следует отметить работы, которые посвящены решению контактных задач о скольжении периодических систем из жестких инденторов с различной формой по вязкоупругому материалу [12], которые позволяют описать поверхности с искусственно заданным макрорельефом (текстурой), а также могут служить моделями для изучения основных закономерностей гистерезисного трения. Трехмерные контактные задачи, рассматривающие скольжение/качение одного индентора по вязкоупругой подложке, описываемой различными моделями материала, в основном решаются с использованием функции Грина и метода граничных элементов [6, 13––18]. Дискретный контакт твердых тел с двоякопериодическим рельефом и вязкоупругого основания рассматривался в работах [19]; контакт жестких волнистых поверхностей с вязкоупругим телом исследован в работах [20, 21]. В двумерной постановке [25] решение получается аналитически с использованием подхода, аналогичного предложенному в известной статье Хантера [22]. Подробный обзор работ по задачам контактного взаимодействия вязкоупругих тел при наличии шероховатости, включающий работы как по индентированию, так и по скольжению, предложен в [23].

Все приведенные выше работы относятся к контакту вязкоупругих тел без покрытия. Во многих случаях рассматривается вязкоупругий слой. Решение для единичного контакта вязкоупругого полупространства (пространственная модель материала), скрепленного со слоем, обладающим изгибной жесткостью, получено в работе [24].

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

2. Постановка задачи и метод решения. Рассматривается скольжение с постоянной скоростью V жесткой периодической системы неровностей под действием среднего давления p_n по поверхности вязкоупругого полупространства с покрытием толщиной H. Покрытие моделируется слоем, обладающим изгибной жесткостью. Используется подвижная декартова система координат, центр которой связан с точкой первоначального касания произвольно выбранной неровности с поверхностью покрытия.

Рис. 1. Модель периодической поверхности n = 2 (a), n = 4 (b), A/L = 1/3.

Скольжение происходит вдоль оси 0x. Форма периодической системы неровностей (рис. 1) описывается следующей функцией:

 $f\left(x,y\right)=\left\{\begin{array}{l}A-A\cdot \left|\cos \left(\frac{x\cdot \pi }{L}\right)\right|\cdot \left|\cos \left(\frac{y\cdot \pi }{L}\right)\right|,n=\mathrm{1,}\\ A-A\left(1-\sin {\left(\frac{x\cdot \pi }{L}\right)}^n\right)\cdot \left(1-\sin {\left(\frac{y\cdot \pi }{L}\right)}^n\right),n=\mathrm{2,}\mathrm{4,}\end{array}\right.$ (2.1)

где A – высота неровностей, L – период неровностей.

Слой и полупространство жестко сцеплены z = – H:

 $w^{\left(1\right)}=w^{\left(2\right)},u_x^{\left(1\right)}=u_x^{\left(2\right)},u_y^{\left(1\right)}=u_y^{\left(2\right)}.$ (2.2)

Здесь w⁽ⁱ⁾, u_x⁽ⁱ⁾, u_y⁽ⁱ⁾ – нормальные и касательные перемещения слоя (i = 1) и полупространства (i = 2).

Свойства линейно-вязкоупругого материала описываются оператором Вольтерра:

 $\left[F\tau \right](t)=1+{\displaystyle \underset{-\infty }{\overset{t}{\int }}\tau (t)}K\left(t-\xi \right)d\xi ,K(t)=k{e}^{(-t/\eta )},$ (2.3)

где t – напряжения, 1/k – время релаксации, h – время последействия.

На верхней границе двуслойного полупространства выполняются следующие условия контакта z = 0:

 $\begin{array}{c}{w}^{\left(1\right)}\left(x,y\right)=f\left(x,y\right)-D,\left(x,y\right)\in {\omega }_i,i=\mathrm{1..}\infty ,\\ {\sigma }_z^{\left(1\right)}=\mathrm{0,}\left(x,y\right)\notin {\omega }_i,i=\mathrm{1..}\infty ,\\ {\tau }_{rz}^{\left(1\right)}={\tau }_{\theta z}^{\left(1\right)}=\mathrm{0,}\forall \left(x,y\right).\end{array}$ (2.4)

Здесь D – осадка неровностей, w_i – область контакта i-й неровности, s_z, s_rz, t_qz – нормальные и касательные напряжения.

Нагрузка на каждую неровность и номинальное давление p_n связаны соотношением:

 $Q_p=p_n{L}^2.$ (2.5)

Уравнение равновесия рассматривается в следующем виде:

 $Q_p={\displaystyle \underset{{\omega }_i}{\iint }p(x,y)dxdy,i=\mathrm{1..}\infty .}$ (2.6)

Для приближенного решения поставленной задачи используется упрощенная постановка:

 $\begin{array}{c}{w}^{\left(1\right)}\left(x,y\right)=f\left(x,y\right)-D,\left|x\right|,\left|y\right|<3L/2\wedge \left(x,y\right)\in {\omega }_i,i=\mathrm{1..9,}\\ {\sigma }_z^{\left(1\right)}=\mathrm{0,}\left|x\right|,\left|y\right|<3L/2\wedge \left(x,y\right)\notin {\omega }_i,i=\mathrm{1..9,}\\ {\sigma }_z^{\left(1\right)}=p_n,\left|x\right|,\left|y\right|\ge 3L/\mathrm{2,}\\ {\tau }_{rz}^{\left(1\right)}={\tau }_{\theta z}^{\left(1\right)}=\mathrm{0,}\forall \left(x,y\right).\end{array}$ (2.7)

В сущности, переход к данной упрощенной постановке означает использование метода локализации [25]. Справедливость подобного перехода доказана в работе [25] применительно к упругому полупространству. Очевидно, что в общем случае для произвольных значений параметров, описывающих вязкоупругое основание, справедливость подобного перехода доказать нельзя. Для каждого конкретного сочетания входных параметров задачи проверка погрешности перехода проводилась в процессе расчетов: если нагрузка на центральную неровность, определяемая при решении контактной задачи, получалась равной Q_p, замена периодической задачи на упрощенную постановку (2.7) полагалась правомерной. Осадка D при переходе от задачи с граничными условиями (2.4) к приближенной задачи с условиями (2.7) становится техническим результатом расчета и при анализе результатов не рассматривается.

Задача с условиями (2.7) решается методом граничных элементов. Для этого используется решение задачи о равномерно распределенном давления внутри квадрата размером 2a Ч 2a, движущегося с постоянной скоростью V по поверхности двуслойного вязкоупругого полупространства [24]. В случае упругого материала вертикальные смещения границы слоя, обладающего изгибной жесткостью, определяются следующим соотношением:

 $w^{\left(1\right)}\left(x,y\right)=-\frac{a}{G}{\displaystyle \underset{0}{\overset{\pi /2}{\int }}{\displaystyle \underset{0}{\overset{\infty }{\int }}\Delta \left(\gamma ,\phi ,\lambda \right)}}\cos \left(\left(x/a\right)\gamma \cos \phi \right)\cos \left(\left(y/a\right)\gamma \sin \phi \right)d\gamma d\phi .$  (2.8)

Здесь g, j – координаты в пространстве двойного преобразования Фурье, l = H/a – безразмерная толщина слоя, G – модуль сдвига. Функция $\mathrm{\Delta }\left(\gamma ,\phi ,\lambda \right)$ определяется из уравнений, полученных в результате двойного преобразования Фурье граничных условий на поверхности покрытия, а также на границе покрытия и полупространства. При этом общее решение представляется в форме Галеркина [26]. Функция $\Delta \left(\gamma ,\phi ,\lambda \right)$ имеет громоздкий вид и линейно зависит от давления в пространстве двойного преобразования Фурье:

 $\overline{q}=q\frac{4}{{\pi }^2}\frac{\sin (\gamma \cos \phi )\sin (\gamma \sin \phi )}{{\gamma }^2\sin \phi \cos \phi }.$ (2.9)

Решение для вязкоупругого полупространства и покрытия получим в результате применения оператора Вольтерра (2.3) к упругому решению (2.8) подобно тому, как это сделано в работе [17].

 $\begin{array}{c}{w}^{\left(1\right)}(x,y\mathrm{,0})=-\frac{a}{2G}{\displaystyle \underset{0}{\overset{\pi /2}{\int }}{\displaystyle \underset{0}{\overset{\infty }{\int }}\Delta (\gamma ,\phi ,\lambda )\cos ((y/a)\gamma \sin \phi )}}\times \\ \times \left(\cos ((x/a)\gamma \cos \phi )+{\displaystyle \underset{-\infty }{\overset{0}{\int }}K(-\tau )\cos (((x/a)+(V/a)\tau )\gamma \cos \phi )d\tau }\right)d\gamma d\phi .\end{array}$ (2.10)

В случае экспоненциального ядра ползучести (2.3) выражение (2.10) примет вид:

 $\begin{array}{c}w(x,y\mathrm{,0})=-\frac{a}{2G}{\displaystyle \underset{0}{\overset{\pi /2}{\int }}{\displaystyle \underset{0}{\overset{\infty }{\int }}\Delta (\gamma ,\phi ,\lambda )\cos \left(\left(y/a\right)\gamma \sin \phi \right)}}\left(\cos \left(\left(y/a\right)\gamma \cos \phi \right)+\begin{array}{c}\\ \end{array}\right.\\ +\left.c\frac{\eta (V/a)\gamma \cos \phi \sin \left((x/a)\gamma \cos \phi \right)+\cos \left((x/a)\gamma \cos \phi \right)}{1+{\left((V/a)\eta \gamma \cos \phi \right)}^2}\right)d\gamma d\phi ,\\ c=k\cdot \eta .\end{array}$  (2.11)

Соотношение (2.11) было использовано для оценки погрешности замены граничных условий (2.4) на условия (2.7). Для этого оценивалась разница вертикальных смещений, возникающих в результате действия постоянной нагрузки вне области | x |, | y | < 3L/2, в центре выделенного квадрата со смещениями от действия сил Q_p, расположенных в узлах квадратичной решетки. Использование факта малости смещений относительно периода решетки позволило показать, что погрешность пропорциональна комплексу $p_{n}c/(G\lambda \eta V),$ то есть чем больше скорость скольжения, тем меньше погрешность замены (2.4) на (2.7). Для решения задачи в постановке (2.7) квадратная область | x |, | y | < 3L/2 разбивается на N элементов, давление в каждом элементе принимаются постоянным. Внутри выбранной области с помощью (2.11) рассчитывается вертикальные смещения w_i⁽¹⁾ границы двуслойного основания под действием номинального давления –p_n, затем с учетом этих смещений, граничных условий, а также условия равновесия для девяти неровностей составляется СЛАУ:

 $\left[\begin{array}{cccc}s& \cdots & s& 0\\ k_1^1& \cdots & k_N^1& -1\\ ⋮& \ddots & ⋮& ⋮\\ k_1^N& \cdots & k_N^N& -1\end{array}\right]\cdot \left[\begin{array}{c}{p}_1\\ ⋮\\ p_N\\ D\end{array}\right]=\left[\begin{array}{c}Q\\ f_1\\ ⋮\\ f_N\end{array}\right]+\left[\begin{array}{c}0\\ {\overline{w}}_1\\ ⋮\\ {\overline{w}}_N\end{array}\right],$ (2.12)

где k_i^j – внедрение в центре элемента j под действием единичного давления в элементе i, которое рассчитывается с использованием формулы (2.11). Решение системы дает распределение давления такое, что граница двуслойного основания повторяет форму неровностей, при этом некоторые давления имеют отрицательные значения, что противоречит условию контакта (2.7), поэтому отрицательные давления приравниваются нулю и система решается заново. Описанная итерационная процедура завершается в том случае, если в решении не окажется отрицательных давлений.

3. Анализ результатов. В данном разделе приведен анализ результатов расчетов, полученных для вязкоупругого полупространства без покрытия и с покрытием – слоем, обладающим изгибной жесткостью. Для анализа использовалась следующая система безразмерных параметров:

 $V^{\prime }=\frac{V\cdot \omega }{L},c=\eta \cdot k,Q^{\prime }=\frac{Q}{E\cdot L^2},A^{\prime }=A/L,H^{\prime }=H/L,\left(x^{\prime },y^{\prime },z^{\prime }\right)=\frac{\left(x,y,z\right)}{L}.$

В результате расчетов определялось безразмерное распределение контактного давления p¢(x, y) = p(x, y)/E, где Е – длительный модуль упругости, и коэффициент трения, обусловленного гистерезисными потерями в вязкоупругом материале. Поскольку в общем случае корректность использования метода локализации для задачи о скольжении периодической системы по вязкоупругому полупространству не доказана, для каждого расчета проводилась численная оценка. Она заключалась в следующем: на центральную неровность должна действовать нагрузка, равная Q_p, с погрешностью не более 10%.

Рис. 2. Распределение контактного давления под девятью неровностями при n = 1, c = 10 A′ = 1/3, (на девять неровностей), V ′ = 8.20 (a), V ′ = 82.05 (b), V ′ = 328.22 (c).

Рис. 3. Распределение контактного давления под девятью неровностями при n = 4, c = 10, A′ = 1/3, Q ′ = 1.23 (на 9 неровностей), V ′ = 8.20 (a), V ′ = 82.05 (b), V ′ = 328.22 (c).

Результаты, представленные на рис. 2–5, получены для вязкоупругого полупространства без покрытия. На рис. 2 и 3 приведены распределения контактного давления, действующего на вязкоупругое полупространство под девятью неровностями, для разных скоростей и различной формы поверхности контр-тела, за которую отвечает параметр n в (2.1). Для используемых параметров расчетов чем больше скорость, тем меньше разница в распределении давлений в системе из девяти неровностей. Это обусловлено более упругим характером контактирования при относительно больших скоростях скольжения (рис. 2, с и 3, с). При увеличении скорости происходит также уменьшение суммарной площади контакта. На неровности, являющиеся передними по направлению скольжения, действуют нагрузки, превышающие среднюю на неровность.

Таблица 1. Распределение нагрузки по неровностям, относительно Q_p, при различных скоростях, n = 1

В табл. 1 показано, как распределяется нагрузка между неровностями (номера неровностей проставлены на рис. 2, а) для различных скоростей скольжения. Входные параметры расчета соответствуют рис. 2. Уменьшение скорости приводит к увеличению погрешности расчетов.

При n = 4 (рис. 3) неровности являются более плоскими, в связи с чем влияние реологических свойств на распределение контактного давления под каждой из неровностей является визуально более выраженным. При этом распределение нагрузки между неровностями (табл. 2) примерно соответствует результатам, полученным при n = 1.

Таблица 2. Распределение нагрузки по неровностям, относительно Q_p, при различных скоростях, n = 4

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

Рис. 4. Распределение контактного давления при (a), n = 4, c = 10, A′ = 1/3, V ′ = 8.20 (a, b), V ′ = 383.22 (c, d), Q ′ = 1.23 (на девять неровностей), Q ′ = 0.137 (на 1 неровность); изолированная неровность (черные кривые), множественный контакт (красные кривые).

На рис. 4 представлены распределения контактного давления (сечение по y ¢ = 0), полученные для центральной неровности в условиях множественного контакта, и для изолированной неровности, на которую действует нагрузка Q_p. Во всех случаях ширина площадки контакта сужается за счет влияния соседних неровностей, что находится в полном соответствии с результатами для упругого полупространства [25]. Для относительно низких скоростей распределение давления существенно несимметрично, особенно для более плоских неровностей, что приводит к появлению момента касательной силы (деформационной составляющей силы трения). При высоких скоростях эффект взаимного влияния увеличивает максимальное значение контактного давления, что соответствует контакту, близкому к упругому. При относительно низких скоростях имеем противоположный эффект за счет большей вытянутости площадки контакта в направлении, перпендикулярном направлению скольжения.

Рис. 5. Зависимость коэффициента трения от скорости для изолированной неровности (a), множественного контакта (б): c = 10, A′ = 1/3, Q ′ = 1.23 (на 9 неровностей), Q ′ = 0.137 (на одну неровность), n = 1 (красная), n = 2 (зеленая), n = 4 (синяя); (c) – сравнение результатов для единичного контакта (сплошная кривая) и множественного контакта (штриховая кривая) при n = 2.

Были получены зависимости, отражающие влияние скорости на деформационную составляющую силы трения при разных значениях параметра n (рис. 5). Следует отметить, что в этих результатах нет монотонности, коэффициент трения для n = 2 больше, чем для n = 4 и 1. Это выполняется как для случая множественного контакта, так и для изолированной неровности. Эффект влияния соседних неровностей можно оценить по рис. 5, с. При скольжении изолированной неровности потери энергии, связанные с несовершенной упругостью, больше. Ранее было показано [17, 18], что зависимость коэффициента трения от скорости является немонотонной – сначала увеличивается, затем уменьшается. В данном случае существенная погрешность замены граничных условий (2.4) на условия (2.7) при околонулевых скоростях не позволила провести соответствующие расчеты.

Рис. 6. Зависимость коэффициента трения от скорости для множественного контакта (a): H ′ = 3.3 · 10^–6, c = 10, A′ = 1/3, Q ′ = 1.23 (на девять неровностей), Q ′ = 0.137 (на одну неровность), n = 1 (красная), n = 2 (зеленая), (синяя); (b) – сравнение результатов для единичного контакта (сплошная кривая) и множественного контакта (штриховая кривая) при n = 2.

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

Рис. 7. Распределение контактного давления при n = 1 (a), n = 4 (b, с): H ′ = 3.3 · 10^–4 (a, b), c = 10, A′ = 1/3, V ′ = 8.20, (на девять неровностей), Q ′ = 0.137 (на одну неровность); изолированная неровность (черные кривые), множественный контакт (красные кривые); (с) – распределения контактного давление при: n = 4 (красная кривая), H ′ = 3.3 · 10^–3 (черная кривая).

Для случая n = 2 проведено сравнение зависимостей коэффициента трения от скорости, полученных для изолированной неровности и для случая множественного контакта (рис. 6, b). Эффект взаимного влияния является более существенным при наличии жесткого изгибающегося слоя, при том что абсолютные значения коэффициента трения при прочих равных условиях больше для вязкоупругого полупространства без покрытия. Влияние соседних неровностей на распределение контактного давления (рис. 7, a, b) при относительно малых скоростях в целом такое же, как и для случая вязкоупругого полупространства, хотя степень не симметрии несколько меньше. Увеличение толщины слоя приводит к увеличению максимального значения контактного давления. При этом следует отметить, что увеличение толщины слоя в десять раз не меняет радикально размер площадки контакта и вид функции p¢(x, y), видимо, в связи с тем, что слой при этом остается относительно тонким.

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

Работа выполнена при финансовой поддержке РНФ (грант 23-19-00484).