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

УДК 539.3+678.065.001.2

А. Е. Белкин, О. А. Одинцов

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

Приведены формулировки и результаты решения контактных задач обжатия неподвижной автомобильной шины и медленного качения обжатой шины по твердой опорной поверхности. Задачи решены методом конечных элементов. Для выполнения контактных условий применен метод штрафа, реализуемый путем модификации матрицы жесткости и вектора узловых нагрузок. Рассмотрен общий случай существования областей сцепления и скольжения. Показано, что в случае скольжения матрица контактной жесткости теряет свойство симметрии. В качестве примера приведены результаты расчетов легковой радиальной шины 175/70R13 Бело-церковского шинного завода.

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

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

Простейшие модели, основанные на представлении шины в виде кольца на упругом основании [1], на сегодняшний день используются редко ввиду недостаточной информативности. Более информативными являются оболочечные модели, например модели оболочки Тимошенко [2], трехслойной оболочки [3, 4], многослойной оболочки [5]. Другая группа моделей основана на соотношениях трехмерной теории упругости. Эти модели позволяют получить информацию о распределении полей напряжений и деформаций в объеме шины. Однако корректная реализация трехмерных моделей предъявляет зачастую невыполнимые требования к вычислительным ресурсам.

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

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

Лагранжа, метода штрафа, либо расширенного метода Лагранжа [6]. Каждый из этих методов в общем случае имеет свои достоинства и недостатки. Метод множителей Лагранжа требует введения дополнительных переменных, однако обеспечивает квадратичную сходимость к точному выполнению условий контакта; метод штрафа обеспечивает квадратичную сходимость без введения дополнительных переменных, однако полученное решение не является точным; расширенный метод Лагранжа объединяет достоинства предыдущих двух методов ценой потери квадратичной сходимости.

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

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

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

Решение позволяет получить распределение контактного давления в пятне контакта, распределение касательных сил, поля перемещений точек каркаса и брекера, НДС слоев и нагрузочные характеристики, полную вертикальную нагрузку Qz, продольную силу Qy и вращающий момент Мх.

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

Предположим, что в пределах возможной площади контакта основание — это пологая, гладкая поверхность, симметричная относительно плоскостей хОг, уОг. Уравнение этой поверхности г = f(x,y).

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

г = Кх,у) + Sz.

При расчете шины смещение основания служит параметром нагруже-ния.

Рассмотрим произвольную точ-Рис. 1. Контакт шины с опорной по- ку наружной поверхности протек-верхностью тора шины, имеющую до деформа-

ции радиус-вектор г0 с координатами х0,у0,г0 и перемещение V с компонентами Ух,Уу. Пространственные координаты этой точки протектора шины должны удовлетворять условию непроникания внутрь основания:

г0 + Vz > -(х0 + Ух, у0 + %) + ^. Введя функцию зазоров между контактирующими поверхностями

П = г0 + Vz - -(х0 + Ух, у0 + Уу) - 8z,

можно переписать условие непроникания как требование неотрицательности зазоров:

П > 0.

Функцию, по знаку противоположную функции зазоров,

д = -п = -(х0 + Ух, у0 + Уу) + 8z - г0 - Vz

будем называть внедрением. Очевидно, что эта функция не может быть положительной:

g

= 0, если точка протектора находится в контакте; < 0, если контакт отсутствует.

Для описания условий контакта может быть использована как функция зазоров, так и функция внедрения. Мы будем использовать последнюю.

Линеаризуем функцию внедрения по отношению к малым перемещениям точек протектора:

g = g0 - - Ф°хvx - ф°чих),

(1)

где д = -(х0, у0) + 8:, - г — начальное внедрение, которое имело

д- д-

бы место в случае проницаемой оболочки, фХ = ^Г" о, Ф = ^Т"

дх у ду

— тангенсы углов наклона касательных к опорной поверхности.

Для пологой поверхности величины ф(°, Ф0 в пределах контакта являются малыми, поэтому будем пренебрегать их квадратами и произведениями по сравнению с единицей. С этой точностью определим единичный вектор внешней нормали к опорной поверхности как

п = —ф°л — ф0] +к

и представим функцию внедрения (1) в виде

д = д° - (п,у) (2)

где у = ух • г + Уу • ] + у2 •к.

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

У = и + АН . (3)

Рассматривая протектор как изотропный винклеровский слой, полагаем (

АН = ^ п, (4)

кр

где дп ^ 0 — нормальное контактное давление; кр — жесткость протектора.

Из соотношений (2)-(4) следует связь между контактным давлением и функцией внедрения:

д = д — (п,и) — т-.

кр

Контактное давление должно быть найдено так, чтобы внедрение обратилось в ноль. Таким образом,

= \ К (д0 - (п и)), если д0 - (п и) >0; (5)

п [О, если д0 — (п, и) < 0.

Область контакта Qc характеризуется неравенством д0 — (п, и) > 0.

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

Пто^ = П + 1 Ц кр (д0 — (п, и))2 ¿П. (6)

пс

Исходя из принципа стационарности функционала (6), получаем вариационное уравнение

önmod = + J I kp

а

(n, и) (n, 6и) — g0 (n, 6u)

dÜ = 0, (7)

которое в матричной форме может быть записано так:

snmod = гп^ Sin}1 k

{n}{n}T in} - g0in}

dü = 0. (8)

Для решения вариационного уравнения (8) применим метод конечных элементов. Процедура МКЭ, основанная на представлении поля перемещений с помощью вектора узловых перемещений {/} в виде

{и} = [Щ{/} ,

3x1 3хппх1

где [N1 — интерполяционные функции; п — количество степеней свободы дискретизированной задачи [7], приводит уравнение (8) к следующему виду:

/ + [Кс]{/}-{^с} = 0 . (9)

дП (дП дП дП 1 Здесь —— = < -—, ——, ... , —— ^ — вектор, компонентами кото-

д {/} 1д/1 д/2 д/п)

рого являются производные функции полной потенциальной энергии по компонентам вектора узловых перемещений;

[К] = Л kp[N]т{п}{п}т(10)

пс

— матрица жесткости контакта,

{^с} = ]т крд0{п}<т (11)

пс

— вектор дополнительной внешней нагрузки, возникающей при контакте.

Для линейно-упругой задачи

п = 1 {/}т [К ]{/}-{/}т {^ },

поэтому уравнение МКЭ (9) принимает вид

[К] + [КС^ {/} = {^} + (12)

где [К] — матрица жесткости системы до вступления в контакт, {^}

— вектор внешних сил.

В случае геометрически нелинейной задачи в уравнение (12) добавляется слагаемое {^ }, учитывающее нелинейный характер зависимости внутренних сил от перемещений:

([К] + [Кс]) {/} + = {^} + (13)

Как видно из уравнений (12), (13) и выражений (10), (11), используемый способ выполнения контактного условия соответствует методу штрафа [8], причем коэффициентом штрафа служит жесткость протекторного слоя.

Рассматриваемая задача заключает в себе два источника нелинейности: во-первых, геометрическую нелинейность оболочки; во-вторых, нелинейность условий контакта (5). Решение находится по двухуровневой итерационной схеме [9]. Согласно этой схеме, внешний итерационный цикл применяется для пересмотра области контакта и выполнения условия (5). Во внутреннем цикле, предназначенном для численного решения системы нелинейных алгебраических уравнений (13) методом Ньютона, область контакта считается неизменной. Итерационный процесс внутреннего цикла продолжается до достижения необходимой степени малости нормы неуравновешенных узловых сил. Внешний итерационный процесс завершается в том случае, когда с заданной точностью достигается совпадение области контакта между двумя соседними итерациями.

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

Я = Ц <1г^ ъ Л qndQ = Ц кр (д0 — {п}т[N11/}) 6П. (14)

ПП с П По

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

Далее приведены результаты расчета статического нагружения легковой шины ВС-43 Белоцерковского шинного завода (Украина) (рис. 2). Для оценки результатов МКЭ использовалось сравнение с решением задачи в дифференциальной линеаризованной постановке с помощью вычислительного комплекса "Каскад" [3].

Симметрия задачи относительно экваториальной плоскости шины позволила выполнить расчет для половины шины. Для решения задачи использовали сетку размером 120 х 46 элементов в окружном/меридиональном направлениях соответственно со сгущением узлов по окружному направлению вблизи пятна контакта. Рассматривали обжатие шины на плоскость либо на беговой барабан диаметром 1,83 метра. Внутреннее давление в шине принималось равным 0,2 МПа.

На рис. 3 представлена нагрузочная характеристика шины при обжатии до 40 мм. Кривая 2 соответствует случаю обжатия на плоскость, кривая 3 — на беговой барабан. Из графика видно, что наличие ненулевой кривизны контактной поверхности приводит к некоторому снижению статической жесткости шины по сравнению со случаем обжатия

Рис. 2. Схематичное изображение профиля шины BC-43

и

Н

6000

5000

то

3000

2000

1000

У

; /г

"з

10 15 20 25 30 35 мм Прогиб шины

Рис. 3. Нагрузочная характеристика шины:

обжатие на плоскость (1,2) и беговой барабан (3)

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

На рис.4 показано распределение контактного давления при обжатии шины на плоскость. На горизонтальной плоскости отложены координаты (Б,ф) где Б — дуговая координата по меридиану (мм), отсчитываемая от точки борта; ф — окружная угловая координата (градусы).

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

Рис. 4. Распределение контактного давления при обжатии на плоскость силой =3700 Н

Рис.5. Сопоставление результатов решения задачи обжатия. Эпюра деформаций корда в каркасе (а) и брекере (б) (масштабы эпюр различны):

сплошные линии — решение геометрически нелинейной задачи МКЭ; штриховые — решение задачи в дифференциальной линеаризованной постановке

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

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

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

С учетом продольных касательных сил я полный вектор контактной нагрузки д принимает вид

д = я*,™ + я/,

где вектор д касательной к недеформируемой опорной поверхности образован как [ п,{ ] и направлен в сторону движения опорной поверхности.

Обозначим радиусы-векторы текущих точек протектора и брекера как Ту и Ги соответственно; угловую координату точек брекера — Также введем угловую скорость П стационарного качения; скорость движения основания V.

Для определения касательных сил в месте контакта используем концепцию "сцепление-скольжение" [1], [11], согласно которой площадка контакта делится на зону сцепления и зону скольжения (рис. 6).

Тангенциальная составляющая контактной нагрузки определяется по модели сухого (кулонова) трения, представляемой неравенством

Ы ^ ¡Яп , (15)

где ¡л — коэффициент трения продольного скольжения шашки протектора по опорной поверхности.

Деформации сдвига элементов протектора и касательные напряжения возникают из-за различия скоростей ГУ, Ги, с которыми движутся в пятне контакта точки протектора и соответствующие им точки бреке-

Рис. 6. Схема определения сдвигов элемента протектора при стационарном качении

ра. Скорость изменения деформации сдвига протектора 7Р составляет

(Г* - Ти) • Г

ъ =—н-'

где НР — высота протектора.

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

дт др

Тогда

г* - ( Пг?1 + Пдр) • Г

Q dYp = -дф hp

где г — расстояние от точки брекера до оси вращения в недеформиро-ванном состоянии; — вектор продольной касательной к недеформи-рованной поверхности брекера.

Из-за отсутствия взаимодействия между шашками протектора примем начальный сдвиг (сдвиг в точках входа в контакт) равным нулю. Кроме того, в дальнейшем выводе будем предполагать, что высота протектора НР не является функцией окружной координаты.

Тогда текущее значение угла сдвига в точке с окружной координатой р будет равно

f

7р(ф) = hp/

fo

f ^ ди\

Q rv - Г1 + дф)

tdy>, (16)

где ф° — угол входа в контакт.

В области сцепления скорость протектора rV • t совпадает со скоростью движения основания V, поэтому, в приближении пологой опорной поверхности, выражение (16) может быть записано в виде

f

Yp4hp = Дк(ф - фо) - r J ti •td(p - (ut |f -ut |fo) = 2° - ut,

fo

где Як = V/Q — радиус качения шины; ut = u • t — проекция перемещения на тангенциальное направление опорной поверхности; д° = Дс(ф - фо) - r(sin ф - sin фо) + и |f0; в случае, когда размеры области контакта малы по сравнению с радиусом шины,

2° = (R - г)(ф - фо) + ut |fo.

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

1

2

Псц = 1 / hp 0(%ц )2dQ,

Шсц _ JJ ^ (и - д0 дП.

Псц

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

ск _ №

1'Р _ .

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

5Аск _ — • £ (П _ — мп5щ ^ .

Пек Пек

Добавляя энергию сдвига протектора к функционалу потенциальной энергии системы (7), получаем следующее вариационное уравнение:

5П*то(1 _ 5П+ кр(ип — д0)5ипд^+

+ II (Щ — д0)5^ 41 ^кр(ип — д°)5щ(Ю. _ 0.

Псц Пск

Это уравнение может быть записано в матричной форме:

дП + ([Кс] + [Ксц] + [Кск]) {/} — (Ш + {ЗД + {Щ) _ 0 , (17)

д {f} где

[Ксц] = JJ G [N]T{t}{t}T[N]dfi;

Псц

{FU = Л G9* [N]T{t}dQ;

Псц

[Кск] = Jl »kp[NY {t}{n}1 [N ]dQ;

nC

[Щ = уу ^Рд0[Щт[г^а

Следует заметить, что матрица контактной жесткости в области скольжения [Кск] не является симметричной, что создает дополнительные трудности в численной реализации.

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

Из решения задачи определяются продольная сила Qy и крутящий момент Мх, которые нужно приложить к оси колеса, чтобы уравновесить силы трения в месте контакта:

Qy = у у Qy dtt ~ у у Qtdü ,

c c

Mx = (qz • y - Qv • (z - Rd))dQ w / I (qn • y - qt • z)dQ + QyRd,

где Ял — динамический радиус шины (расстояние от оси колеса до опорной поверхности).

Значение вертикальной силы Qz, аналогично задаче статического нагружения, вычисляется по формуле (14).

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

На рис. 7 приведено распределение касательных сил в пятне контакта шины, в котором можно наблюдать области сцепления и скольжения. Для расчета был использован коэффициент трения ц = 0,7 и сетка конечных элементов размером 80 х 46 элементов в окружном/меридиональном направлениях соответственно.

На рис. 8 показаны графики изменения натяжения нитей брекера шины в экваториальной плоскости по окружному направлению вблизи области контакта. График 1 соответствует решению о се симметричной задачи о нагружении шины внутренним давлением 0,2 МПа. График 2 — обжатию неподвижной шины; кривые 3-5 соответствуют различным режимам медленного качения. Для графика 3 область скольжения распространяется на всю область контакта, что соответствует режиму

Рис. 7. Распределение касательных сил в пятне контакта

£cord

О,ООП

0,0010

0,0006

0,0002 О

°'°00-80 -60 -40 -20 0 20 40 60 ср,°

Рис. 8. Влияние касательных сил на деформацию нитей брекера:

1 — нагрузка внутренним давлением; 2 — обжатие неподвижной шины; 3 — режим полного скольжения; 4 — режим полного сцепления; 5 — общий случай медленного качения

буксования колеса. График 4 соответствует режиму полного сцепления шины с основанием, в котором область скольжения отсутствует. График 5 иллюстрирует режим качения, в котором присутствует как область сцепления, так и область скольжения. Из рис. 8 видно, что наличие касательных сил в области контакта и выбор режима качения оказывают влияние на распределение деформаций и напряжений внутри шины.

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

Рис. 9. Схема используемого для расчетов конечного элемента трехслойной оболочки

Конечный элемент трехслойной оболочки. Для моделирования оболочки используется восьмиузловой конечный элемент, схематично изображенный на рис. 9. Вектор узловых степеней свободы размерности 24, включающий в себя линейные перемещения узлов, имеет вид [/} = [«1 VI ... щ у8 ^з}.

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

В пределах каждой из граней используется билинейная аппроксимация перемещений; по толщине прослойки перемещения аппроксимируются линейно, основываясь на аппроксимации перемещений соответствующих точек лицевых граней. Для исключения эффекта "сдвигового заклинивания" прослойки используется схема аппроксимации сдвигов, аналогичная описанной в работе [12].

Матрица тангенциальных жесткостей всего элемента строится из матриц тангенциальных жесткостей мембранных граней и прослойки по принципу суперпозиции:

[К ] =

24x24

■[Kmi] 12x12

0

0

[Km2] i2 i2

+ [Кз] ,

24 24

где [Ктг1] и [Ктг2] обозначают соответствующие матрицы тангенциальных жесткостей мембранных граней; [К3] — матрица жесткостей прослойки.

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

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

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

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

1. Б у х и н Б. Л. Введение в механику пневматических шин. - М.: Химия, 1988, 223 с.

2. Белкин А. Е., Нарская Н. Л. Численный анализ деформаций автомобильной шины при стационарном качении // Сб.: Проблемы прикладной механики, динамики и прочности машин. - М.: Изд-во МГТУ им. Н.Э. Баумана, 2005. - С. 71-89.

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

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

5. Григолюк Э. И., Куликов Г. М., Плотникова С. В. Контактная задача для пневматической шины, взаимодействующей с жестким основанием // Проблемы машиностроения и надежности машин. - 2004. - №4. - C. 55-63.

6. Wriggers P., Wagner W., Stein E. Algorithms for non-linear contact constraints with application to stability problems of rods and shells // Computational Mechanics. - 1987. - №2. - P. 215-230.

7. Зенкевич О. Метод конечных элементов в технике. - М.: Мир, 1975. -541 с.

8. Т у р ч а к Л. И. Основы численных методов. - М.: Наука, 1987. - 318 с.

9. Nour-Omid B., Wriggers P. A two-level iteration method for solution of contact problems // Computer methods in applied mechanics and engineering. -1986. - Vol.54. - P. 131-144.

10. Григолюк Э. И., Шалашилин В. И. Проблемы нелинейного деформирования. - М.: Наука, 1988. - 231 c.

11. Karaoglan L., Noor A. K. Sensitivity analysis of frictional contact response of axisymmetric composite structures // Computers & Structures. - 1995. -Vol. 55. - №6. - P. 937-954.

12. B a t h e K. J., D v o r k i n 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. - P. 697-722.

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