CC BY
12УДК 531/534+681.3.06
Реализация модели Контенсу-Эрисмана касательных сил в контактной задаче Герца
И. И. Косенко, Е. Б. Александров
Российский государственный университет туризма и сервиса 141221, Московская обл., Пушкинский район, Черкизово-1, ул. Главная, д. 99, kosenkoQccas.ru, aleksandrovebQgmail.com
Получено 15 февраля 2009 г.
В рамках контактной задачи Герца строится приближенная модель вычисления результирующего мотора (винта) касательных в контакте сил сухого трения. Винт состоит из суммарной силы трения и момента (пары сил) трения верчения. Рассматриваемый подход естественным образом развивает построенную ранее компьютерную модель герцевого упругого контакта. Силы сухого трения и момент этих сил интегрируются по эллиптическому пятну контакта. В общем случае аналитическое вычисление упомянутых интегралов приводит к громоздким выражениям, составленным из десятков слагаемых, являющихся рациональными функциями, зависящими, в свою очередь, от полных эллиптических интегралов с модулем — эксцентриситетом контактного пятна. Для реализации достаточно быстрой компьютерной модели касательных сил проводится приближенное построение в направлении, предложенном еще Контенсу. Представляемая здесь модель является естественным развитием упрощенной модели Контенсу в следующих направлениях: а) модель является анизотропной — суммарные силы трения вдоль главных осей контактного эллипса в общем случае различны; б) для поступательных и почти поступательных относительных движений в области контакта используется регуляризованный кулоновский закон трения; в) построена также приближенная модель момента трения верчения. Для верификации модели используются результаты, полученные ранее несколькими авторами. В качестве тестового динамического примера используется модель волчка тип-топ. Оказалось, что процесс «переворота» волчка на сферу меньшего радиуса («голову»), численно смоделированный при помощи подхода, основанного на применении техники многозначных отображений, практически совпадает с численной верификацией представляемой здесь модели. Динамическая модель шарикоподшипника используется для детального сравнительного тестирования различных подходов к вычислительной реализации касательных сил. Объекты модели упругих контактов между шариками подшипника и его внутренним и внешним кольцами, основанные на законе Кулона касательных сил точечного контакта, были заменены с учетом описываемого здесь модифицированного подхода Контенсу. Оказывается, упрощенные формулы подхода Контенсу обеспечивают скорость моделирования даже бблыную, чем в модели точечного контакта.
Ключевые слова: контактная модель Герца, модель Контенсу-Эрисмана, упрощенная модель Контенсу, модель В. Г. Вильке, волчок тип-топ, модель шарикоподшипника
1.1. Kosenko, Е. В. Aleksandrov Implementation of the Contensou-Erismann model for tangent forces in the
Hertz contact problem
In frame of the Hertz contact problem an approximate model to compute resulting wrench of the dry friction tangent forces is built up. The wrench consists of the total friction force and the drilling friction torque. An approach under consideration develops in a natural way the contact model built up earlier. The dry friction forces and torque are integrated over the contact elliptic spot. Generally an analytic computation of the integrals mentioned leads to the cumbersome calculation, decades of terms, including rational functions depending in turn on complete elliptic integrals. To implement the elastic bodies contact interaction computer model fast enough one builds up the approximate model in the direction as it was proposed by Contensou. The model under construction is one derived from the Contensou simplified model in the following directions: (a) the model is anisotropic: the total friction forces along ellipse axes are different; (b) for the translatory and almost translatory relative motions one uses the Coulomb friction law regularization; (c) the approximate model for the drilling torque also has been constructed. To verify the model built the results obtained by several authors were used. The Tippe-Top dynamic model is used as a an example under testing. It turned out the top revolution process is identical to one simulated using the set-valued functions approach. The ball bearing dynamic model is used to verify different approaches to the tangent forces computational implementation in details. The model objects corresponding to contacts between the balls and raceways were replaced by ones of new class developed here. Then the old friction model of the regularized Coulomb type and the new one, approximate Contensou, each embedded into the whole bearing dynamic model were thoroughly tested and compared. It turned out the simplified Contensou approach provides the computer model even faster in compare with the case of the point contact.
Keywords: the Hertz contact model, the Contensou-Erismann model, the Contensou simplified model, the Vilke model, the Tippe-Top, the ball bearing model
Mathematical Subject Classification: 00A71 68U20 T0E18 T0E55 T4H15 74M10 74M15
1. Постановка задачи
Решение задачи Герца [1] для нормального давления, распределенного по площадке контакта эллиптической формы, задается [2] при помощи формулы
а (x,y) =
3N
1 _ ^
2irab V а2 Ь2 ’
где Ж — величина суммарной силы нормального давления, а Ь — большая и малая полуоси наружного эллипса контактной площадки (рис. 1.1), Рхуг — локальная система координат контакта, ориентированная так, что ось х направлена вдоль большой полуоси эллипса.
Рис. 1.1. Площадка области контакта.
На рис. 1.1 предполагается, что тело А располагается «снизу» (иод чертежом), а тело В — «сверху», и все вычисляемые здесь силы действуют на тело В со стороны тела А. Рассмотрим методику вычисления силового винта, составленного из Е = ¥х\ + — касатель-
ной силы и Т = Тгк — крутящего момента результирующей пары, возникающих вследствие действия сил сухого трения, распределенных по поверхности контактного эллипса.
По определению сухого трения, касательное напряжение т(х,у) в точке М(х,у) контактной площадки вычисляется по формуле
т (x,y) = -fa (x,y)
V м(х,у) |vM(x,y)| ’
(1.1)
где f — величина коэффициента трения, vm(x,y) — скорость относительного проскальзывания тела B то телу A в геометрической точке M. Формула (1.1) используется при vm(x, y) = 0. При vm(x, y) = 0 вектоp т(x, y) не определен — он может иметь произвольное направление и иметь длину не более величины fa(x,y). По этой причине для описания задач с ударами и сухим трением часто применяются формализмы многозначных (set-valued) отображений и функций [3]. Попробуем построить вычислительные процедуры для задач с сухим трением, оставаясь в рамках классического анализа и опираясь на известные результаты [4] об асимптотической близости точной и приближенной задач.
AB
ми, то возможны следующие варианты относительного движения в контакте.
В первом случае угловая скорость ^ = uzк относительного верчения равна нулю. От-
B
скоростей относительного проскальзывания по контактной площадке, задаваемое известной формулой Эйлера
7 vpx - Uz y
VM(x,y)= Vp + VM/p(x,y)= | vpy + Uzx I , (1.2)
где (рис. 1.1) Ум/р(ж, у) = [ш, г (ж,у)\, г (ж, у) = ж1 + yj — радиус-вектор текущей ТОЧКИ М (ж, у) контактной площадки, является постоянным и равным у р.
Чтобы вычислить главный вектор Е сил сухого трения, следует в соответствии с моделью Контенсу-Эрпсмана [5, 6] найти интеграл по эллиптической площадке контакта
Е = // т (ж,у)«у. (1.3,
В нашем случае этот интеграл сводится к виду
2паЬ |ур|УУ V а2 62 |ур|'
Здесь снова при Ур = 0 вектор Е не определен, принимая любые возможные значения из круга радиуса fN.
Е
хого трения
|ур |
ранее успешно применявшуюся [7] в случае точечного контакта твердых тел. Здесь функция д$ (у) задается по формуле
0 при у ^ 0,
д5 (у) = { у/6 при 0 < у ^ 5,
1 при у > 5.
Величина 5, характеризующая переход от режима относительного скольжения к режиму относительного почти качения и обратно, выбирается достаточно малой, но такой, чтобы обеспечить сохранение надежности вычислительного процесса.
Главный момент сил сухого трения относительно центра приведения Р имеет вид
Т = //[г(*.у),т (ж,у)\«у. (1.4,
Поскольку векторы г (ж, у) и т (ж, у) находятся в плоскости Ржу, то вектор Т имеет только одну ненулевую компоненту — проекцию на ось Ря. Очевидно, эту проекцию можно вычислить по формуле
Т = JJ (жту(ж, у) - утх(ж, у)) ^Ыу.
Это выражение при = 0 имеет вид
3/^ 1 ГГ,............... [ ж2 у2
2шь Ы 11^ \А - ¥ - Ь'Ыу-
из которого в силу свойств антисимметрии при отображениях ж ^ -^1 у ^ —у следует, что
Т = 0.
Рассмотрим теперь случай, когда угловая скорость относительного верчения отлична от нуля. В этом случае Ум (ж, у) = 0 всюду в плоскости Ржу, кроме одной точки — мгновенного центра скоростей, который может находиться либо внутри, либо вне контактной
площадки. В любом случае знание векторной функции т(ж, у) в этой точке никак не влияет на результат интегрирования по формулам (1.3), (1.4), поскольку точка имеет нулевую меру Лебега, а функция т(ж, у) в ней ограничена величиной /а(ж, у). Поэтому всюду в
ЕТ вать выражение (1.1), считая, что Ум (ж, у) = 0 и = 0.
2. Предельные случаи
Представляется интересным прежде всего рассмотреть два предельных случая вычисления интегралов (1.3), (1.4), соответствующих часто встречающимся примерам упругого контактирования в инженерных приложениях.
Замечание 2.1. Справедливости ради следует заметить, что интегралы (1.3), (1-4) вычисляются аналитически в конечном виде. Однако получающиеся выражения содержат десятки слагаемых, зависящих в виде рациональных функций от параметров задачи и от полных и неполных эллиптических интегралов первого и второго рода, в свою очередь зависящих от этих параметров. В таком громоздком виде задача представляется малопригодной для компьютерной реализации, и тем более — в приложениях реального времени, при моделировании контактного взаимодействия в системах, состоящих из достаточно большого количества тел.
Возвращаясь к упомянутым выше предельным случаям, нетрудно убедиться [8], что при большой динамической жесткости контактирования (большой модуль упругости Юнга материалов тел) полуоси контактного эллипса являются величинами, малыми по сравнению с радиусами главных кривизн поверхности относительного контакта, получаемой вычитанием ¿-координат поверхностей тел в локальной системе координат. В этом предельном случае малой является величина параметра а ^ 1.
Кроме того, в большом количестве инженерных приложений в контактных задачах преобладает относительное качение тел, а верчение либо считается почти отсутствующим
~ 0), либо является очень незначительным ^ 1). Такая ситуация реализуется, например, в подшипниках качения.
Параметры а, ^ неудобны тем, что они являются размерными величинами и зависят от выбора единиц измерения. Перейдем к более удобному безразмерному малому параметру
Ур Р
где у — величина относительной скорости скольжения в точке Р, а ф — угол, составленный Рж Ур
ным переменным интегрирования (£, п) по формулам
где е — эксцентриситет контактного эллипса. Преобразование (ж, у) ^ (£,п) задаваемое
1
Так что уравнение (1.2) удобно представить в форме
ур = у = а = соя ф, в = йш ф,
(2.1)
Ум(£, п) = у + е
а — £7п
в + 0
где величина
аш
V
и есть искомый безразмерный параметр, малый в двух упомянутых выше предельных случаях. Далее по примеру работы [9] вместо размерной будем использовать безразмерную
скорость относительного скольжения в точке Р: и = £-1.
Используя новые переменные интегрирования, получим новые выражения для компонент главного вектора и главного момента сил трения
Г, = с,[[ , ^ (22)
У У д/1 — 2е{аг)Г) — /?£) + е2 (£2 + 72т/2)
^ II (2-3)
Т = с [[ У1 - Р - ^ [~(а74 - + * «2 + (2.41
г УУ \ 2$(а~/г) — /?0 + в2 '- + 72|/2) ’
где размерные множители вычисляются по формулам
3/Ж 3/Жа
2п 2п
Здесь снова, как и выше в случае скорости (но без изменения обозначений), вместо размерных будем использовать безразмерные величины сил и моментов, нормируя их при помощи размерных значений /Ж, /Жа соответственно. Тогда вместо размерных множителей с/, с^ далее будем использовать выражения
3 3
С/ 2тг ’ С* 2тг'
Нетрудно заметить, что подынтегральные функции в (2.2-2.4) представимы в виде степенных рядов (равномерно сходящихся в пределах некоторого круга конечного радиуса в С) по £. С учетом условия £ ^ 1 ограничимся для определенности и простоты первыми слагаемыми в этих рядах. Прежде всего, рассмотрим подынтегральный множитель
11
у/1 - 2є(агут) — /?£) + є2 (£2 + 72т/2) л/1 — 2ssmv + s2'
где приняты обозначения
грі—¡—ñ~ñ . OL'yrj — /3£
s = £д/ч + 7Т) smz/
причем, как легко заметить, V — угол между вектором vp и РМ.
Получилось хорошо известное выражение для производящей функции полиномов Лежандра [10]. Так что
л/1 — 2ssinz/ + s2
J2snpn (sin V),
n=0
причем, первые три полинома от sin v имеют вид
3 1
Ро (sin v) = 1, Pi (sin v) = sin v, P2 (sin v) = - sin2 г/ — -.
Отсюда получаем требуемое разложение
1
д/1 — 2e{a^r¡ — /?£) + е2 (£2 + 72т/2)
= 1 + е(а7п - в£)+
+e2
+
Введем полезные для дальнейшего обозначения
1ш,п = 7" ff Ví-e- v4mVndCd7]
Тогда из (2.2)-(2.4) имеем
Fo
F
F
i
3 r
w/0)0,
3
“7Г“ [(j + 2/?w)/i>0 - (i + 2/3w)/o i] ,
2П
2
To
Ti
3
2vr
l2,o - (3a^w - aj - ei)1i,i +
+
3 2 1 ,
-a;--------w — ш
2 2
¿0,2
— 7Г- (Ph,o — &Io,i) k,
3
— — [(l — (32) /2,0 + 2a/3Iiti + (l — «2) /0,2] k.
(2.5)
(2.6)
(2.7)
(2.8)
(2.9)
(2.10)
(2.11)
Остается вычислить интегралы вида (2.5). Чтобы привести эти интегралы к стандартному виду, в качестве новых переменных интегрирования можно использовать полярные координаты: £ = р cos p, n = Р sin Р- Тогда
i
2п
1 = Yn 1m,n — /
J д/1 — р2 pm+n+l dp I i У cos™ p sin” pdp
,0 / \0
Используя свойства четности функций, нетрудно убедиться, что
2п
cosm p sinn pdp = 0,
o
если хотя бы одно из целых неотрицательных чисел то, п нечетно. Так что в формулах (2.72.11) отличны от нуля только величины
2п 2п 2п72
-*0,0 = -*2,0 = ту, -*0,2 = - •
3 15 15
Поэтому
Fo = -w,
Fi = 0,
F
T
T
1
2 =-----<----W — в"
5 12
-а w — ш
2 2
3. Упрощенная модель
Задаем направление вектора vp так, что угол ф будет фиксирован. Проведем исследование компонент силы Е и момента Т как функций от скалярной переменной и ^ 0 величины вектора относительной скорости в центре контактной площадки. В соответствии с (2.2)-(2.4) имеем точные выражения
2п
Fx(u) Fy(u)
ó г--2 j a-ejpsmp
oo
i 2n
— — J д/1 — P2pdp J -
/3 + ep eos p
IC(p,¥>)l
dp,
(3.1)
(3.2)
3
1
0
i
i
i
3
Tz(u) = —— V1 - p2p2dp x
.........2Л2,
2vr
o
2n
f cos p(e + ep cos p) - Y sin p(a - eYp sin p)
x J------------------kmi--------------------dv' m
o
где Z(p, p) = VM (p cos p,YP sin p) и в подынтегральных выражениях знаменатель дроби имеет вид
|С(р,р)| = \Jl — 2ep(ajsinр — f3cos p) + e2p2 (cos2 p + y2 sin2 p).
Нетрудно заметить, что разложение функций от u, задаваемых при помощи (3.1-3.3), e
точки u = го по степеням 1/u. Удерживая в разложениях (2.6) члены наинизшего порядка, имеем
F°°(u) = Ф§° + ^Ф^°, Т°°(п) = “Т^°, (3.4)
где следует положить Ф0 = Fo, Ф0 = F2, TO = Ti. Сравнительные графики функций Fx(u) и FX°(u), Fy(u) и Fj0(u), Tz(u) и T0(u) показаны на рис. 3.1, 3.2, 3.3 соответственно. Здесь для определенности положено: ф = п/4, e = 0,81. Графики точных функций отображаются сплошной линией, а их аппроксимаций — пунктиром.
Рис. 3.1. Асимптотика ж-координаты вектора силы при и ^ ж.
Рис. 3.2. Асимптотика у-координаты вектора силы при и ^ ж.
Рис. 3.3. Асимптотика момента трения верчения при и ^ ж.
Найденная при и ^ ж асимптотика представления силы и момента сухого трения дает достаточно хорошие приближения для точных величин искомых интегралов в широком
ии (аш в размерных единицах) теперь уже не малый параметр е стремится к единице, и разложение для производящей функции полиномов Лежандра становится расходящимся. Аппроксимация при помощи представлений (3.4) уже не будет иметь хорошего качества. Не говоря о том, что при и ^ 0 функции го (3.4) не ограничены Е^(и) ^ ж, Тте(и) ^ ж, в то время как точные интегралы в (3.1)-(3.3) дают при и = 0 конечные (нулевые в (3.1)-(3.2)) значения, в чем можно убедиться непосредственным вычислением.
Чтобы выправить ситуацию, рассмотрим на этот раз асимптотику интегралов при и ^ 0. Вычисления показывают (для проверки можно использовать [11, с. 172-178]), что
ГХ(и) = — аАи + 0(и2), Гу (и) = — вВи + 0(и2), Тх (и) = —С + 0(и2), (3.5)
причем после вычисления в (3.1-3.3) повторных определенных интегралов (сначала по р, а затем по р) получим
в = |(к(е) + ВД^ОД), с = |Е(е),
где К(е) Е(е) — полные эллиптические интегралы первого и второго рода. Заметим, что при комбинировании рассматриваемого здесь алгоритма вычисления винта сил сухого трения с ранее реализованным алгоритмом вычисления нормальной силы в контактной задаче Герца [8] перечисленные выше полные эллиптические интегралы уже будут вычислены в процессе расчета нормальной силы. Поэтому здесь, при вычислении винта касательных сил, для этого не потребуется никаких дополнительных вычислительных ресурсов.
Как было замечено еще Контенсу [5], главным эффектом в модели сухого трения Контенсу-Эрисмана при наличии контактной площадки, наряду с возникновением крутящего момента пары сил трения верчения, является монотонное убывание до нуля величины главного вектора сил трения.
С другой стороны, из (3.5) легко заметить, что при аш ^ 0 «крутизна» изменения результирующей силы (как функции от ь) растет. При этом либо контактная площадка стремится к точке, либо исчезает угловая скорость верчения. В конце концов, при достаточно малой величине аш, когда А, В ^ аш/£, мы будем иметь уже реализованный ранее [7] для точечного контакта случай регуляризации модели сухого трения. Таким образом, результирующая компьютерная модель всегда будет использовать «регулярный» случай силы трения, убывающей, хотя и достаточно «круто» к нулю.
и
от кулоновской на величину второго порядка малости, вслед за Контенсу [5] упростим нашу модель, считая силу трения кулоновской при u € Щ,у, го) и линейной при u € [0>и£,у). Заметим, что в «поведении» силы здесь имеет место анизотропия: постоянные и£ и иу вдоль осей Рж и Ру, вообще говоря, отличаются. Найдем эти величины из уравнений
АиХ = 1, Виу = 1.
Теперь приближенную модель Контенсу для функций РХ(и), Гу(и) Тх (и) можно представить в виде
„ , ч Г Аи при и € [0, иХ),
хо (и) = — а< * ч
[ 1 при и € [иХ, го),
уг (и) = _в / Ви ПРИ и € [°>иу)>
С ( ) в \ 1 при и € [иу, ГО ,
Т (и) = [ — С при и € [0,и*) ,
С() \ Т~(и) при и € [и*, го),
где и* является точкой склейки горизонтальной «полки» высоты С и ветви убывающей на бесконечности гиперболы, задаваемой функцией Т^(и). Уравнение, задающее значение величины и*, имеет вид
с =
и%
Сравнение графиков функций точной (сплошная) и построенной здесь приближенной (пунктир) моделей Контенсу-Эрисмана представлено на рис. 3.4, 3.5, 3.6 при значениях параметров, использованных для построения рис. 3.1, 3.2, 3.3.
Рис. 3.4. ж-координата вектора силы точной и приближенной моделей.
Рис. 3.5. у-координата вектора силы точной и приближенной моделей.
Рис. 3.6. Момент трения верчения точной и приближенной моделей.
Если все найденные величины п£, п*, и* в процессе движения становятся меньше 5/аи, то мы приходим к использовавшейся ранее регуляризованной модели Кулона. Так что построенная здесь приближенная модель Контенсу-Эрисмана является простейшим обобщением упомянутой кулоновской, еще более ее регуляризуя добавлением параметров п£, и*, и*, расширяющих область применимости параметра регуляризации 5. Такая регуляризация есть просто следствие наличия контактной площадки в точной модели. Ясно, что она не
требует (за исключением случаев а = 0 ш = 0) использования многозначных отображений, позволяя нам применять технику классического анализа.
Построенная в данной работе приближенная модель имеет несколько существенных отличий от кусочно-линейной аппроксимации, построенной в статье [5]:
а) модель анизотропна и пригодна для эллиптических контактных площадок любого эксцентриситета;
б) при мгновенном поступательном и почти поступательном относительном движении тел, когда одновременно выполняются условия
и%аш < 5, и*аш < 5, и*аш < 5, применяется регуляризация сухого трения, предложенная в [4] и в [12];
в) для момента трения верчения также используется приближенная модель.
Вычисления показывают, что представленные здесь выражения дают более точное приближение модели Контенсу-Эрисмана, чем дробно-линейная аппроксимация, удовлетворяющая граничным условиям в нуле и на бесконечности. Применение в аппроксимациях Паде полиномов второй и третьей степени [13] улучшает качество приближения, но требует уже значительно большего объема вычислений.
4. Случай модели В. Г. Вильке
Как известно [8], хорошее качество аппроксимации нормальной упругой силы контактного взаимодействия в достаточно широком диапазоне эксцентриситетов дает формула В. Г. Вильке. Компьютерная модель такой реализации контактного взаимодействия работает заметно быстрее реализации точной модели Герца. Главной причиной этого является необходимость в последнем случае в процессе моделирования решать трансцендентное уравнение вида
О(с) = д, 0 < с < 1, 0 <д < 1 (4.1)
относительно с = е2 — квадрата эксцентриситета контактной площадки, функция левой части задается равенством
°{с)=1Щ)~{1~с)> (4-2)
где в качестве аргумента полного эллиптического интеграла первого рода, как и в [14], для
с
Алгоритм В. Г. Вильке для вычисления нормальной контактной силы не требует знания с
бы не терять полученный при расчете нормальной силы выигрыш и не тратить, теперь уже при расчете касательной силы трения, компьютерное время на решение уравнения (4.1), упростим это уравнение так, чтобы довести процесс его решения до однократного применения явной линейной формулы.
Для этого исследуем более подробно, чем в [8], поведение левой части уравнения (4.1) — функции О (с). Известно [8], что функция С (с) при с € [0,1] монотонно убывает от единицы до нуля. Причем О'(с) ^ —ж при с ^ 1. Попытаемся приблизить О (с) на отрезке [0,1] конечным числом ломаных. Оказалось, что это можно выполнить при помощи двух сегментов
с
вый, качество аппроксимации, «накрывает» область эксцентриситетов, близких к единице, когда сам алгоритм В. Г. Вильке уже не обеспечивает хорошего качества приближения нормальной упругой силы. Так что в «рабочем» диапазоне эксцентриситета этого алгоритма нам все же удается с высокой точностью заменить О (с) линейной функцией от с.
Пусть с1 € (0,1) является параметром, характеризующим кусочно-линейную аппроксимацию О* (с) функции О(с) так, что при с € [0, с1 ] аппроксимирующая функция вычисляется по формуле
С* (с) = 1 — —с,
где коэффициент линейного члена выбирается из очевидного условия
-!=
Второй линейный участок функции О* (с) протяжен от точки с = с* до точки с = 1. Так что на нем аппроксимация имеет вид
С* (с) — — — с)-
1 — с1
Качество приближения функции О(с) (сплошная линия) при помощи О* (с) (пунктир) при с1 = 0,4 можно оценить на рис. 4.1.
Рис. 4.1. Функция О (с) и ее аппроксимация.
В результате применения алгоритма вычисления нормальной упругой силы в контакте [8] в каждый текущий момент моделирования рассчитывается локальный базис контакта и коэффициенты Р, Q канонической формы квадратичной части представления разности уравнений поверхностей взаимодействующих тел в локальной системе координат в виде
г = Рж2 + Qy2, Р < Q.
Полагая С* = Р/^, находим величину с* — приближенное значение квадрата эксцентриситета эллиптической площадки контакта из уравнения
О* (с) = С*.
После вычисления с* для расчета упомянутых выше величин А В С находим полные эллиптические интегралы первого и второго рода при помощи тэта-функций [14]. Прежде всего, при любом с* € [0,1) можно использовать быстро сходящееся разложение
ГО
0з(д) = 1 + 2£; д”2
”=1
так, что полный эллиптический интеграл первого рода будет вычисляться по формуле
К (с) =
где параметр д вычисляется с очень высокой точностью при помощи соотношений [14, 15]
_ 1 1 - (1 - с)1/4
£ “ 2 ' 1 + (1 - с)1/4’
д = е + 2е5 + 15е9 + 150е13 + 1707е17 + ..., причем выписанных в разложении величины д членов достаточно, чтобы обеспечить точ-
1/221
Для полного эллиптического интеграла второго рода удобно использовать формулу [16]
Е(с) = -К(с) + ^
2”
со
--2У —--------2
12 ¿К1- ч2п)
К(с) с
Е(с) — К (с) _ <Ж 2с ^с
для регуляризации выражения левой части можно использовать хорошо в этом случае сходящиеся гипергеометрические разложения [17].
5. Вычислительные эксперименты
Верификация представленной здесь модели касательных сил проводилась в два этапа: а) для случая кругового контакта; б) для случая некругового эллиптического контакта. В первом случае рассмотрена динамическая модель известного волчка тип-топ. Все параметры, определяющие эту механическую систему, совпадают с данными, представленными для этого примера в статье [18], авторы которой, в свою очередь, взяли эти параметры из более ранней работы [19].
Геометрически «твердое» тело волчка составляется из двух шаров (рис. 5.1): большего, радиусом К = 1,5 ■ 10-5 м, и меньшего, радиусом г = 0,5 ■ 10-5 м. Центр масс всего волчка считается расположенным «под» центром большего шара на оси симметрии на расстоянии ао = 3 ■ 10-3 м и на расстоянии а1 = 16 ■ 10-3 м от центра меньшего шара. Масса волчка равна т = 6 ■ 10-3 кг. Главные центральные моменты инерции таковы (тело предполагается динамически симметричным): экваториальный момент инерции составляет 8 ■ 10-7 кг - кг2, а полярный — 7■ 10-7 кг -кг2. Материал, из которого изготовлен волчок, и горизонтальный пол, по которому он катится, предполагается деревом с модулем упругости Юнга Е = 9,1 х
\ в
[ г\ Ъ
\аЛ
1 Л \
\\«° \ К
Рис. 5.1. Геометрические параметры волчка.
х109 Н/м2, что соответствует (при значениях коэффнцнентов Пуассона в 0,3) эффективному модулю упругости Е* = 5 • 109 Н/м2. Коэффициент трения сколвжения предполагается равнвш / = 0,3.
В начальный момент времени центр масс волчка считается неподвижным, а сам волчок (точнее — его больший шар) никак не вдавлен в упругую горизонтальную поверхность, меньший шар расположен в верхней полусфере шара большего, а начальный угол наклона осп симметрии волчка составляет 9о = 0,1 рад. Начальная угловая скорость также в точности совпадает с данными работы [18], равна ^0 = 180 с-1 и направлена вдоль осп динамической симметрии.
Заметим, что в работе [18] контактная задача интерпретируется в ограниченном смысле: площадка контакта считается заданной, и в этом предположении вычисляются вектор силы трения и момент трения верчения, а именно: круговой контакт вычисляется для случая статического равновесия упругой силы со стороны плоскости и силы тяжести волчка. На самом же деле в процессе движения сила нормальной реакции, реализуемой здесь при помощи силы упругости, меняется. Следовательно, меняется и площадка (радиус) пятна контакта. Так что волчок должен испытывать также вертикальные упругие микровибрации. При этом, как можно видеть из изложенного выше в упрощенной модели, будет меняться момент трения верчения, что, в общем случае, может отразиться на движении волчка.
В модели, построенной в данной работе, рассматривается неограниченная контактная задача: упругая нормальная сила вычисляется из модели Герца (либо В. Г. Вильке) с добавлением (нелинейного) вязкого слагаемого. Одновременно вычисляется эллиптическое пятно контакта, и эти данные используются для расчета касательной силы и момента трения верчения в рамках упрощенной модели Контенсу.
Замечательным образом, вычислительный эксперимент показал, что сценарий перевертывания волчка с «ног» (опора на больший шар) на «голову» (опора на меньший шар), полученный в [18] при помощи совершенно иной вычислительной модели, с высокой точностью воспроизводится и у нас. График рис. 5.2 хорошо иллюстрирует этот процесс: так же, как и в [18], на нем представлена эволюция угла наклона 9 (£). Сравнение этого графика с соответствующим участком эволюции угла наклона в [18] подтверждает их полную идентичность. Единственное отличие состоит в наличии на рис. 5.2 дополнительного «дрожания», что свидетельствует о высокочастотной составляющей в динамике упругого контакта волчка и плоскости. Аналогичную идентичность показывают и кривые индикаторов контактов шаров, из которых составлен волчок, с плоскостью (рис. 5.3). Индикатор контакта большого шара представлен сплошной линией, а малого — штрихованной линией.
Рис. 5.2. Эволюция оси симметрии волчка.
0.024
0.0 0.8 1.6 2.4 3.2 4.0
Рис. 5.3. Эволюция индикаторов контакта. Интервал переворота увеличен.
Помимо описанной была проведена также верификация при помощи результатов работы [9]: формулы для силы трения и момента трения верчения, полученные в [9] (случай а) были применены в компьютерной модели динамики волчка, построенной на языке Modélica в варианте неограниченной контактной задачи. В рамках той же компьютерной модели была реализована описываемая здесь упрощенная модель Контенсу (случай б) и модель, использующая дробно-линейную аппроксимацию Паде для сил и момента трения (случай в). Результаты эволюции угла наклона показаны на рис. 5.4. При этом случаю а соответствует сплошная кривая, случаю б— штрихованная линия, случаю в — пунктир. Легко видеть, что сценарии переворота волчка ближе всего друг к другу в случаях ап б. Снимок трехмерной анимации волчка тип-топ показан на рис. 5.5.
3.5
1----'----1---'----1---'----1---'----
0 12 3 4
Рис. 5.4. Сравнительный анализ трех моделей.
Рис. 5.5. Снимок трехмерной анимации волчка.
Замечание 5.2. Поведение, характерное для волчка тип-топ (переворот «па голову») можно легко получить и в рамках модели «регуляризоваппого» трения Кулона. Регуляризацию при этом следует понимать в смысле, предложенном в работах |4, 12| и использовавшемся в случае точечного контакта в |7|. «наклоняя» график зависимости силы трения от скорости в окрестности пуля так. чтобы заменить этим наклоном разрыв первого рода в пуле. Чем более пологим является наклон, тем скорее обнаруживается динамический эффект переворота волчка. Как показывает приближенная модель Коптепсу, именно этот наклон появляется в соответствующем графике зависимости силы трения от скорости в рамках уже точной модели Коптепсу Эрисмапа.
На втором этапе верификации, в случае существенно эллиптической контактной площадки. рассматривалась, как и в работе [8]. динамическая модель шарикоподшипника. При этом сравнивались два подхода: а стандартная модель нормальной силы Герца плюс упрощенная модель сил трения Контенсу; б упрощенная модель нормальной силы В. Г. Вильке плюс упрощенная модель сил трения Контенсу. Как и в работе [8]. динамические модели случаев а и б мало отличаются друг от друга. При этом модель б демонстрирует большую (на 20 %) производительность. В качестве иллюстрации на рис. 5.6 сравнивается одна из составляющих касательной силы трения в контакте между одним из шариков подшипника и его внутренним кольцом. Случаю а соответствует сплошная кривая, а случаю б штри-
Рис. 5.6. Эволюция одной из компонент силы трения. Увеличен интервал завершения процесса моделирования.
хованная. Вычислявшиеся в процессе моделирования значения квадрата эксцентриситета контактной площадки в этом примере оказались постоянными по времени и равными 0,687 в случае а и 0,643 — в случае б.
6. Обсуждение результатов
Подводя итоги вышеизложенному, отметим следующие основные выводы.
• Эффект переворота волчка тип-топ «па голову» полностью обусловлен «регуляризацией» силы сухого трения, обеспечиваемой моделью Контенсу-Эрисмана, в окрестности нулевой скорости проскальзывания. Вычислительные эксперименты показали, что если в модели Кулона наклон «регуляризованной» силы трения в окрестности нуля достаточно велик, то эффект «тип-топ» либо не наблюдается совсем, либо, кратковременно реализуясь в течение долгого времени эволюции движения, быстро исчезает. И только заметное уменьшение упомянутого наклона немедленно приводит к перевороту волчка «на голову» с дальнейшей длительной прецессией в этом положении.
Именно такое свойство силы трения и дает модель Контенсу-Эрисмана. При этом роль момента трения верчения сводится к диссипативному эффекту последующего постепенного «падения» волчка, приближающего его к статически устойчивой конфигурации.
•
Контенсу, фактически уже найдены в рамках алгоритма Герца для расчета нормальной силы, то с вычислительной точки зрения применение этой модели практически «ничего не стоит». Если же еще учесть эффект регуляризации, доставляемый моделью Контенсу-Эрисмана, то мы приходим к неожиданному, на первый взгляд, результату: численное моделирование динамики контактного взаимодействия при помощи мо-
делей Герца (нормальная сила) и Контенсу-Эрисмана (касательная сила и момент трения верчения) оказывается на 20-30 % быстрее, чем комбинирование моделей Герца и «простого» сухого трения Кулона. Замедление счета при использовании модели связано, конечно же, с большой жесткостью задачи в этом последнем случае.
• Хотя для изотропного случая контактирования средние значения касательных сил моделей Кулона и Контенсу-Эрисмана не слишком сильно различаются, однако в анизотропном случае первая модель становится неадекватной, в то время как вторая продолжает обеспечивать «правильное» моделирование процесса контактирования. Это свойство становится важным, например, в случае симуляции шариковых подшипников, в которых контактные площадки имеют заметно вытянутую эллиптическую форму.
Один из авторов выражает благодарность П. С. Красильникову за полезные обсуждения. Работа выполнена при поддержке грантов РФФИ 08-01-00600-а, 08-01-00718-а, 08-08-00553-а.
Список литературы
[1] Hertz H. Uber die Berührung fester elastischer Körper // J. Reine Angew. Math., 1882, vol. 92, pp. 156-171.
[2] Ландау JI. Д., Лифшиц E. M. Теоретическая физика: В 10 тт.: Т. 7: Теория упругости. М.: Наука, 1987. 248 с.
[3] Leine R.I., Nijmeijer H. Dynamics and bifurcations of non-smooth mechanical systems. Berlin-Heidelberg-New York: Springer, 2004. 236 pp.
III Новожилов И. В. Условия застоя в системах с кулоновым трением // Изв. АН СССР. МТТ, 1973, № 1, с. 8-14.
[5] Contensou P. Couplage entre frottement de glissement et frottement de pivotement dans la théorie de la toupie // Kreiselprobleme Gydrodynamics: IUTAM Symp. Celerina. Berlin: Springer, 1963, pp. 201—216. [Контенсу П. Связь между трением скольжения и трением верчения и ее учет в теории волчка // Проблемы гироскопии. М.: Мир, 1967, с. 60-77.
[6] Erismann Th. Theorie und Anwendungen des echten Kugelgetriebes //Z. Angew. Math. Phys., 1954, vol. 5, pp. 355-388.
[7] Косенко H. H. Реализация компьютерной модели динамики систем твердых тел с освобождающими связями // Матем. моделирование, 2006, т. 18, № 2, с. 95-106.
[8] Александров Е. Б., Вильке В. Г., Косенко И. И. Контактная задача Герца: Численная редукция и объемометрическая модификация // Журн. вычисл. матем. и матем. физ., 2008, т. 48, № 12, с. 2195-2211.
[9] Журавлев В.Ф. О модели сухого трения в задаче качения твердых тел // ПММ, 1998, т. 62, вып. 5, с. 762—767.
[10] Бейтмен Г., Эрдейи А. Высшие трансцендентные функции: Функции Бесселя, функции параболического цилиндра, ортогональные многочлены. М.: Наука, 1974. 296 с.
[11] Градштейн И. С., Рыжик И. М. Таблицы интегралов, сумм, рядов и произведений. М.: Физмат-гиз, 1962. 1100 с.
[12] Rooney G. T., Deravi P. Coulomb friction in mechanism sliding joints // Mech. Mach. Theory, 1982, vol. 17, no. 3, pp. 207-211.
[13] Kireenkov A. A. Three-dimensional model of combined dry friction and its application in non-holonomic mechanics //Proceedings of ENOC-2005, 5th EUROMECH Nonlinear Dynamics Conference (Eindhoven, The Netherlands, August 7—12, 2005) / Eds: D.H.van Campen, M. D. Lazurko, W. P. J. M. van den Oever. Eindhoven: Eindhoven Univ. of Technology, 2005, pp. 571577.
[14] Уиттекер, Э. Т., Ватсон, Д. H., Курс современного анализа. Ч. П. Трансцендентные функции Москва: Физматгиз, 1963, 516 с.
[15] Янке E., Эмде Ф., Лёш Ф. Специальные функции. М.: Наука, 1964. 344 с.
[16] Milne-Thomson L. М. Elliptic integrals // Handbook of mathematical functions: with formulas, graphs, and mathematical tables / Eds: M. Abramowitz, I.A.Stegun. New York: Dover, 1972, pp. 589-626.
[17] Бейтмен Г., Эрдейи А. Высшие трансцендентные функции: Эллиптические и автоморфные функции, функции Ламе и Матье. М.: Наука, 1967. 300 с.
[18] Leine R.I., Glocker Ch. A set-valued force law for spatial Coulomb-Contensou friction // Eur. J. Mech. A Solids, 2003, vol. 22, no. 2, pp. 193-216.
[19] Friedl C. Der Stehaufkreisel: Master’s thesis. Univ. Augsburg, Inst, für Physik, 1997. 60 p.
