Конечно–элементный анализ контакта автомобильной шины с опорной поверхностью на основе оболочечной модели Текст научной статьи по специальности «Физика»

где А^ — износ пары трения за одно включение на г-м режиме работы. Допустимый износ дисков трения определяется сопоставлением ограничений, налагаемых конструкцией фрикциона и самого диска. Для МКД, работающих всухую, он лимитируется толщиной металло-керамических накладок и способом их крепления к основе диска.

СПИСОК ЛИТЕРАТУРЫ

1. Жучков В. М., Филиппов А. Н., Юдин Е. Г. К расчету гидродинамических потерь мощности в дисковых фрикционных узлах // Известия вузов. Машиностроение. - 1982. - № 4. - С. 3-7.

2. Крагельский И. В. Трение и износ. - М.: Машиностроение, 1968. - 480 с.

3. Трение, изнашивание и смазка: Справочник / Под ред. И.В. Крагельского и В.В. Алисина. - М.: Машиностроение, 1979. - Т. 2. - 358 с.

4. Дроздов Ю. Н. Обобщенные характеристики для оценки износостойкости твердых тел // Трение и износ. - 1980. - Т. 1. - № 3. - С. 417-424.

5. Фанталов В. С., Жучков В. М., Юдин Е. Г. Износостойкость фрикционных элементов, работающих в масле // Известия вузов. Машиностроение. -1982. - № 12. - С. 59-62.

Статья поступила в редакцию 24.02.04

УДК 539.3

А. Е. Белкин, Н. Л. Нарская

КОНЕЧНО - ЭЛЕМЕНТНЫЙ АНАЛИЗ КОНТАКТА АВТОМОБИЛЬНОЙ ШИНЫ С ОПОРНОЙ ПОВЕРХНОСТЬЮ НА ОСНОВЕ ОБОЛОЧЕЧНОЙ МОДЕЛИ

Рассмотрены формулировка и основные этапы решения геометрически нелинейной контактной задачи обжатия шины. Для выполнения условий контакта использован метод конечных элементов с процедурой множителей Лагранжа. Деформации шины описываются теорией оболочек типа Тимошенко. Расчеты выполнены с помощью четырехугольного четырехузлового оболочечного элемента М1ТС4, для которого подробно представлены определяющие соотношения и результаты тестирования. Приведен пример расчета обжатия автомобильной шины.

Исследование контакта автомобильной шины с опорной поверхностью является одной из наиболее важных задач механики шин. Численному решению этой задачи посвящены работы [1-6], в которых используются различные по степени детальности моделирования и точности подходы. В основном используются две расчетные модели: сло-

истой оболочки [1-3] и трехмерного тела [4-6]. Обе модели реализуются, как правило, методом конечных элементов (МКЭ).

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

В настоящей работе шина представлена как упругая слоистая оболочка с конечной жесткостью при поперечном сдвиге. В отличие от работ [2, 3], где контактная задача решалась с помощью алгоритма Уд-завы, далее применен метод множителей Лагранжа, обеспечивающий более быструю сходимость итерационного процесса.

Формулировка и метод решения контактной задачи. Задача о контакте упругого тела (шины) с жестким основанием формулируется [7] как задача минимизации функционала полной потенциальной энергии

П (и) = У а (и) ¿V д ■ и <т (1)

на множестве функций перемещений и, удовлетворяющих в области возможного контакта Ое кинематическому условию непроникания точек тела внутрь основания:

Пп - ип > 0, (2)

где а (и) — функция удельной энергии упругой деформации; V — объем тела; д — внешняя нагрузка на части поверхности тела ; пп — начальные зазоры между деформируемым телом и основанием в направлении нормали п к основанию; ип = и ■ п — нормальные составляющие векторов перемещений.

Искомое контактное давление должно удовлетворять условию неположительности

Чп 6 0. (3)

В области возможного контакта условия (2), (3) могут быть объединены:

Пп - ип > 0 ,

Чп 6 0,

Чп (Пп - ип) = 0, при Я е

здесь Я — радиусы-векторы точек наружной поверхности тела до деформации.

Применение процедуры МКЭ приводит к дискретным аналогам выражений (1), (2). В случае линейно-упругой системы функция полной потенциальной энергии имеет вид

П ({Д}) = 1 {А}т [K] {А} - {А}т {f} ,

где {А} — вектор перемещений в узлах сетки конечных элементов; [K] — матрица жесткости; {f} — вектор узловых сил.

Условие непроникания (2) в i-м узле сетки конечных элементов записывается следующим образом:

{с}т {А} - Пиг 6 0, (4)

где {сг} — известный числовой вектор; пиг — значение зазора для этого узла.

Для решения контактных задач наиболее часто применяют метод множителей Лагранжа и метод штрафов, позволяющие минимизировать функцию потенциальной энергии при наличии ограничений. В настоящей работе использован метод множителей Лагранжа. Для каждого ограничения (4) вводится множитель Аг, имеющий смысл контактной силы в i-м узле, и составляется модифицированная функция потенциальной энергии:

m

П mod ({А} , {А}) = П ({А}) + ^ Аг ({Сг}т {А} - Пиг) =

г=1

= 1 {А}т [K] {А} - {А}т {f} + {А}т ([C]т {А} - {п}), (5)

где {А} = {А! ... Ат} — вектор множителей Лагранжа; {п} = = {Пп1 ... Ппт} — вектор дискретных зазоров; [С] = [{сх} ■ ■ ■ {ст}] — матрица, составленная из векторов, определяющих кинематические условия (4) в узлах; т — общее число кинематических условий.

Условия стационарности этой функции по переменным {А}, {А} приводят к системе уравнений

[К] {А} + [С] {А} = {/} , [С]т {А} = {п} , (6)

решая которую получают перемещения {А} и контактные силы {А}. Далее проверяется условие неположительности контактных сил, и по результатам проверки пересматривается область контакта. Уравнения (6) приходится формулировать и решать многократно для различных пробных зон контакта.

При практической реализации метода для регуляризации задачи к функции Лагранжа (5) добавляют слагаемое (1/2к) {а}т {а} [8]. В результате второе уравнение системы (6), описывающее условия точечного контакта, заменяется следующим:

[с]т {А} + (1/к) {а} = {п} .

Применительно к задаче обжатия шины, решаемой по оболочечной модели, смысл добавленного слагаемого (1/к) {а} очевиден — это деформации протектора.

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

[к ]. ¿{А}+1 + [с ] {аь+1 = {/} - {£ },;

[с]т¿{АЬ+1 + (1/к){аь+1 = {п} - [С]т {А},; {А},+1 = {АЬ + ¿{А},-+1, з = 0, 1, 2, ■■■ ,

где {А},, ¿{А},+1 — приближенные значения перемещений и малые поправки к ним; {#}, — вектор внутренних сил при перемещениях {А},; д {5 }

— матрица тан-

= д {А} генциальной жесткости при перемещениях {А},.

Оболочечная модель ши-

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

Рис. 1. Профиль легковой радиальной шины 175/70К13

Э.И. Григолюка и Г.М. Куликова [9]. На основе всестороннего изучения осесимметричной задачи о надувке шин авторами сделан вывод о том, что модель оболочки с глобальной (не послойной) аппроксимацией перемещений по толщине может применяться в качестве первого приближения, предназначенного для оценки общей нагруженности материалов шины без детального анализа локальных эффектов в конструктивно сложных зонах борта и кромок брекера. В настоящей работе применение теории оболочек типа Тимошенко рассматривается именно в указанном смысле.

Примем срединную поверхность каркаса шины в качестве базовой поверхности оболочки и определим ее деформации {г} = г2 712} соотношениями квадратичной теории:

£1 = ех + ±К2 + 0?); = е2 + 2К2 + 02); 712 = + + 0102, (7)

где е1, 01 (1 ^ 2) — параметры деформации поверхности, линейно зависящие от перемещений.

Деформации эквидистантных слоев {г}г = { £1г г2г 712г } определяются добавлением деформаций изгиба и кручения:

{г}, = {£} + г {к}, (8)

где {к} = { к1 к2 2к12 } — вектор, составленный из производных углов поворота вь в2 материальной нормали к базовой поверхности оболочки [9]. Значения поперечных сдвигов {7} = { 713 723 }, осред-ненные по толщине оболочки, вычисляются как разности углов поворота материальной и геометрической нормалей оболочки:

713 = в1 - 01; 723 = в2 - 02. (9)

Деформационным соотношениям (7)-(9) соответствуют нелинейные уравнения равновесия, которые после приведения напряжений к базовой поверхности и аппроксимации перемещений с использованием узловых переменных {А} получают следующий вид:

[ ([£Г]т{т} + [В]т{м} + [в7]т{д}) ^ = {/}, (10)

где {Т}, {д} — векторы мембранных и поперечных сил; {М} — вектор моментов в оболочке; [В^], [В7], [Вк] — матрицы связи вариаций деформаций оболочки с возможными узловыми перемещениями, т.е. матрицы из соотношений

ад = [В£^]5{А}; 5 {к} = [ВК]5{А}; 5{7} = [В ]5{А}. (11)

В написании матрицы мембранных деформаций [В^] верхние значки напоминают о нелинейном характере зависимости этих деформаций от перемещений.

При формулировке связи между внутренними силами {Т}, {М}, {д} и деформациями {г}, {к}, {7} предполагаем, что и резина, и ре-зинокордные слои проявляют свойства линейной упругости в пределах деформаций, наблюдаемых при нормальных условиях эксплуатации шин. В этом случае физические соотношения для оболочки имеют вид

{Т} = [Я£]{г} + [ЗДк};

{М} = [Я£кМ + [£к]{к}; (12)

М = [Я7 ]{7 };

где [Яе] = У [Е]^г; [Б^] = У г[Е]^г; ^ г2[Е]^г;

н н н

[Д7] = У [С] ¿г; [Е], [С] — матрицы модулей упругости резины и ре-

н

зинокордных слоев при плоском напряженном состоянии и поперечных сдвигах.

Уравнения (7)-(12) составляют принципиальную основу рассматриваемой оболочечной модели.

Получим выражение для матрицы тангенциальной жесткости оболочки. Для этого найдем приращение вектора внутренних сил {д} (т.е. левой части уравнений равновесия (10)), соответствующее малому приращению вектора узловых перемещений:

d{g} = / (d[ßfL]T{T} + [BfL]T4T}+

+ [Бге]т4м} + [В]т¿п. (13)

С учетом соотношений упругости (12) и формул для деформаций (11) выражение (13) приводится к виду

¿{д} = ([к] + [Кк] + [К] + [К]) 4Д} + У {Т(14)

о

где

[K] = [BfL]T[D£ ][BfL]dO;

n

[Кк] = ([BfL]T [d£k ][bK] + [bK]t [d£k ][ßfL0 dO;

[К ] = [Бк]т[^к ][BK]dO;

[K7 ] = [B7 ]т [D7 ][B7 ]dfi

— матрицы мембранной, мембранно-изгибной, изгибной и сдвиговой жесткостей оболочки соответственно.

Исходя из нелинейных формул (7) определения деформаций базовой поверхности, устанавливаем вид матрицы [В^]:

[Ве*£] = [Ве] + [0][В,],

где [Ве] — матрица, определяющая линейную по перемещениям часть мембранных деформаций; [В#] — матрица связи вектора углов поворота {0} = { 01 <^2 02 } элементов базовой поверхности с век-

01 0 0

тором узловых перемещений;

0 0 0 0

матрица

углов поворота.

Теперь последнее слагаемое выражения (14) может быть преобразовано следующим образом:

d[ßfL]x{T = [В^[0]т{Т =

где [T] =

= у [Б^]т[T]d{tf}dfi = iy ß]T[T]ß]dfi|d{A},

— матрица мембранных усилий.

Ti 0 0 0

0 Ti 0 S

0 0 T2 0

0 S 0 T2

Таким образом, получаем окончательные выражения для приращения вектора внутренних сил

¿Ы = ([Ке] + [Кек] + [Кк] + [К7] + [К]) ¿{А} и матрицы тангенциальной жесткости оболочки

[К = [К ] + [К£к ] + [Кк ] + [к7 ] + [К ],

где [Кст] = / [В#]т[Т][В#— так называемая матрица начальных на-п

пряжений.

Рис. 2. Плоский четырехугольный оболо-чечный элемент с 24 степенями свободы

Оболочечный конечный элемент MITC4. Известно [10], что при расчете оболочек с учетом поперечных сдвигов МКЭ в форме метода перемещений может наблюдаться значительное завышение жесткости конечно-элементной модели, называемое "сдвиговым заклиниванием". Многим оболочечным элементам присущ этот дефект, ограничивающий их применимость. В настоящей работе для расчетов использован элемент MITC4 (Mixed Interpolation of Tensorial Components) [11], у которого ситуация "заклинивания" устранена благодаря специальной процедуре аппроксимации поперечных сдвигов.

MITC4 — это четырехугольный четырехузловой элемент с билинейной аппроксимацией текущих координат и перемещений. Элемент может быть как плоским, так и закрученным, если четвертый узел выходит из плоскости первых трех узлов. Однако для расчета шины, вся базовая поверхность которой может быть разделена на конечные элементы линиями кривизны, достаточно использовать плоские элементы. Ограничимся описанием этого варианта.

Введем локальную декартову систему координат: оси x, y расположим в базовой плоскости элемента, а ось z направим по нормали к этой плоскости. В каждом из четырех узлов элемента определяются шесть степеней свободы: три линейных перемещения щ, v, wi по направлениям осей x, y, z и три угла поворота 6xi, 0yi, Qzi относительно тех же осей (рис. 2). Общее число степеней свободы элемента равно 24.

Заметим, что углы 9x, 9y представляют собой использованные в модели оболочки повороты материальной нормали: вх = —в2, = А.

На элементе проводится билинейная аппроксимация текущих коор-динат

{x У} = Ni(f, n){xi Vi}

(15)

i=1

и перемещений

4

{u v w dx 6y} = Ni (f, П){Ui Vi Wi öyi}, (16)

i=1

1

где N (£, п) = 4(1 + 6 £ )(1 + П П) = 1,4) — функции формы; £, п — локальные нормированные координаты, изменяющиеся в пределах от —1 до +1; £ = ±1, п = ±1 — координаты узлов элемента.

Построение матриц жесткости [Ке], [Кек], [Кк], [Кст] осуществляется обычным для изопараметрических элементов образом.

При формировании матрицы сдвиговой жесткости [К7] поступают так. Деформации поперечного сдвига 713, 723 (для плоского конечного элемента эти деформации далее обозначаются 7жг, 7уг) включаются в расчет с помощью аппроксимации, независимой от аппроксимации деформаций изгиба к1, к2, 2к12. Сначала через узловые перемещения выражаются углы сдвига в плоскостях границ элемента в срединных точках А, В, С, Д (см. рис. 2):

Yz = _ l(e>i + ^)

YS

S12

w4 _ w3

1

_ — + )

s43 2

yD = w1 _ w4 1 (й + , )

YVz =--ö + 4)

s14 2

1

Ynz =--TT (0x2 + 0Ж3)

s23 2

yi _ У2

S12

У4 _ У3 s43

У1 _ У4 s14

У2 _ У3 s23

+ ^ + ^У2)

+ ^ (^У4 + ^y3)

+ 2 (öy1 + öy4)

+ 2 (öy2 + Öy3)

Ж1 _ x2 S12

Ж4 - Ж3

s43 X1 X4

s14 x2 _ x3 s23

(17)

где в^ — длины сторон элемента.

Затем, основываясь на значениях сдвигов в указанных четырех точках, проводится аппроксимация законов изменения сдвигов внутри элемента:

Hz =

Yn2

A

S12Ygz

se

s14Y

D nz

1 + n 2

1+1 2

+

+

с

s437gz

se

s23Y

B nz

1 — n

1 — i

(18)

где , вп — длины соответственно £-, п-линий, являющиеся переменными в произвольном четырехугольнике.

Аппроксимация (18) позволяет описать состояние чистого изгиба элемента, что не удается получить непосредственно из полей перемещений (16). Для реализации чистого изгиба достаточно, чтобы в ноль обращались поперечные сдвиги 7^, 7^, , 7ПВ. При этом деформации изгиба

К1 = , К2 =--—, 2К12 =

дж

д y

dy дх

могут принимать произвольные значения. Полученная возможность предсказывать нулевые или малые значения деформаций сдвига, не

s

s

n

n

ограничивая другие деформации, предотвращает эффект "сдвигового заклинивания".

Итак, формулами (17), (18) определены сдвиги в координатных плоскостях n = const, £ = const. Теперь необходимо перейти к углам сдвига Yxz, Yyz. Для этого используем формулы преобразования деформаций при повороте осей:

Ye* = Yxz cos а + Yyz sin а; Yn* = Yx* cos ф + Yy* sin ф, (19)

где а, ф — углы между осями ж, £ и ж, n соответственно. Обратные соотношения имеют вид

Yxz = (те* sin ф - Ynz sin a)/ sin(^ - а); Yyz = (-Ye* cos ф + Ynz cos а)/ sin(ф - а).

(20)

Направляющие косинусы осей £ и п выражаются через элементы матрицы Якоби преобразования координат (15)

дж д у _

|£ д£

д ж ду

. дп дп -

[J ] =

следующим образом:

дж дy

cos а = A- —; sin а = A- —;

e д£ e д£

. дж дy

cos ф = An —; sin ф = An —,

где A = \l (M) 2 + (D , An Ш2 + (— масштабы

длины £ - и n-линий.

Учитывая значения направляющих косинусов, представим формулы прямого и обратного преобразования деформаций (19), (20) в виде

{AeJez AvyVz} = [J]{y*Z Yyz};

(Yxz Yyz} = [J]-1{AeYez AnYnz}.

-ir, „ , (21)

Чтобы получить окончательную аппроксимацию поперечных сдвигов 7жг, 7уг, заменим в соотношении (21) сначала деформации 7^, 7пг выражениями (18), а затем деформации в срединных точках А, В, С, Б их значениями по формулам (17). Учтем также, что при билинейной

аппроксимации координат (15) длины £-, п-линий связаны с масштаба-

ми длин отношениями sg = 2Ag, действий получим

= 2An. В результате названных

i(e)

где

{ Yxz Yyz } = [J]-1 [N7][ß!]{A7}

{A7 }(e) = {W1|0

— подвектор узловых перемещений элемента, связанных с изгибом;

BY ] = t

2 У21 X12 —2 У21 X12 0 0 0 0 0 0

0 0 0 0 0 0 —2 У34 х43 2 У34 х43

2 У41 х14 0 0 0 0 0 0 —2 У41 х14

0 0 0 2 У32 X23 —2 У32 X23 0 0 0

— матрица связи вектора поперечных сдвигов

— „ „,в

{ s12Ygz S43 Ygz S14Y,,z S23 Y^ }

с узловыми перемещениями {А7}(е);

1 1 + n 1 — n 0 0

2 0 0 1 + i 1 — i

— матрица функций формы для поперечных сдвигов.

Таким образом, для элемента сформирована матрица поперечных сдвигов

[в7] = и-1к][в;].

Дальнейший ход формирования матрицы сдвиговой жесткости соответствует традиционной технологии МКЭ.

Пример расчета. В качестве примера рассмотрим результаты расчета легковой радиальной шины 175/70R13, профиль которой показан на рис. 1. Основные параметры конструкции шины таковы: наружный диаметр 582 мм, диаметр посадки на обод 329 мм, ширина профиля по пресс-форме 178 мм, каркас однослойный, жесткость корда в каркасе 290 Н/нить, частота укладки корда 1320 нитей/м, брекер двухслойный, жесткость корда в брекере 39000 Н/нить, частота укладки корда 640 ни-тей/м.

Для оценки точности решений, получаемых с помощью элемента М1ТС4, проводились тестовые расчеты при осесимметричном нагру-жении шины внутренним давлением. Решение МКЭ сравнивалось с результатами численного интегрирования дифференциальных уравнений

s

n

этой одномерной задачи, кото -рые в рассматриваемом случае служили эталоном.

Расчет МКЭ выполнялся для симметричной четверти шины, базовая поверхность которой покрывалась сеткой 47 х 37 узлов.

На рис. 3 показаны графики изменения вдоль профиля шины деформаций поперечного сдвига (а), а также изменения меридиональных (б) и окружных (в) деформаций в точках наруж-ной поверхности корпуса шины при внутреннем давлении 0,2 МПа. Сложный характер распределения деформаций по профилю объясняется тем, что коэффициенты упругости оболочки [Бе ], [Б£к], [Б к] имеют сильно различающиеся значения в трех зонах: в зоне борта, боковой стенки и брекера. На границах этих зон коэффициенты упругости меняются разрывно, что приводит к скачкообразному изменению меридиональных деформаций и резкому изменению других компонент.

Из графиков видно, что результаты расчета МКЭ и числен -ного интегрирования уравнений согласуются качественно и весьма близки по числовым значениям. Таким образом, элемент М1ТС4 обеспечивает приемлемую точность расчета оболочек по модели Тимошенко.

Перейдем к результатам решения контактной задачи. На рис. 4 показано распределение контактного давления по ширине беговой дорожки при нагрузке на шину 3500 Н. Показанное давление является условным, поскольку получено как отношение узловых сил {Л} к брутто

Рис. 3. Сопоставление результатов расчета МКЭ и численного интегрирования уравнений для осесимметричной задачи надувки шины:

сплошные линии — результаты интегрирования; точки — результаты расчета МКЭ

площади контакта, т.е. без учета расчлененности протектора. Реальное же давление на выступы протектора несколько выше условного, его можно определить, разделив условное давление на коэффициент насыщенности рисунка протектора.

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

В рассмотренном примере в процедуре определения узлов контакта и значений контактных сил сходимость наблюдалась на 6-й итерации метода множителей Лагранжа, при этом число внутренних итераций метода Ньютона равнялось 8. По сравнению с решением этой же задачи методом Удзавы время счета сократилось многократно.

Таким образом, апробированы алгоритм и оболочечный конечный элемент М1ТС4, дающие надежные результаты в нелинейной контактной задаче обжатия автомобильной шины.

Рис.4. Распределение контактного давления по ширине беговой дорожки:

1 — в центральном сечении пятна контакта; 2, 3, 4 — в меридиональных сечениях, повернутых по отношению к центральному на углы 5°; 10о; 15° соответственно

Рис. 5. Деформация корда в каркасе (а) и каркасной резины (б) в направлении, поперечном корду

СПИСОК ЛИТЕРАТУРЫ

1. Николаев И. К. Математическая модель диагональной шины и ее численная реализация // Механика пневматических шин. - М.: НИИШП, 1976. - С. 5-36.

2. Чернецов А. А. Решение контактной задачи для пневматической шины с использованием геометрически нелинейной теории оболочек. Дисс.. ..канд. техн. наук. - М.: МГАДИ, 1993. - 138 с.

3. Белкин А. Е. Разработка системы моделей и методов расчета напряженно-деформированного и теплового состояний автомобильных радиальных шин. Дисс.... д-ра техн. наук. - М.: МГТУ им. Н.Э. Баумана, 1998. - 284 с.

4. Rothert H., Idelberger H., J a c o b i W., L a g i n g G. On the finite element solution of the three-dimensional tire contact problem // Nuclear Engineering and Design (North-Holland), 1984. - Vol. 78. - PP. 363-375.

5. F a r i a L.O., Oden J.T., Yavari B., Tworrydlo W.W., Bass J.M., Becker E. B. Tire modeling by finite elements // Tire Science and Technology, 1992. - Vol. 20. - № 1. - PP. 33-56.

6. Ш е ш е н и н С. В., М а р г а р я н С. А. Численное моделирование контактного взаимодействия шины с дорогой // 11-й симпозиум "Проблемы шин и резинокордных композитов". - М.: НИИШП, 2000. - Т. 2. - С. 211-217.

7. Кравчук А. С. Вариационный метод исследования контактного взаимодействия и его реализация на ЭВМ // Расчеты на прочность. Вып. 25. - М.: Машиностроение, 1984. - С. 33-50.

8. Si mo J. C., Wr i g g er s P., T a y l o r R. L. A perturbed Lagrangian formulation for the finite element solution of contact problems // Computer methods in applied mechanics and engineering. 1985. - Vol. 50. - PP. 163-180.

9. Григолюк Э. И., Куликов Г. М. Многослойные армированные оболочки. Расчет пневматических шин. - М.: Машиностроение, 1988. - 288 с.

10. Голованов А. И., Корнишин М. С. Введение в метод конечных элементов статики тонких оболочек. - Казань, Казанский физ.-техн. ин-т, 1989. - 269 с.

11. Bathe K. J., Dvorkin E. N. A formulation of general shell elements — the use of mixed interpolation of tensorial components // International journal for numerical methods in engineering. 1986. - Vol. 22. - PP. 697-722.

Статья поступила в редакцию 12.01.04

Александр Ефимович Белкин родился в 1951 г., окончил в 1974 г. МВТУ им. Н.Э. Баумана. Д-р техн. наук, профессор кафедры "Прикладная механика" МГТУ им. Н.Э. Баумана. Автор более 50 научных работ в области расчетов на прочность, механики пневматических шин.

A.Ye. Belkin (b. 1951) graduated from the Bauman Moscow Higher Technical School in 1974. D. Sc. (Eng.), professor of "Applied Mechanics" department of the Bauman Moscow State Technical University. Author of more than 50 publications in the field of structural analysis, mechanics of pneumatic tires.