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

УДК 539.3+678.065.001.2

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

О.А. Одинцов

МГТУ им. Н.Э. Баумана, Москва, Российская Федерация e-mail: odintsov@bmstu.ru

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

Ключевые слова: метод конечных элементов, гиперупругость, резина, шина.

DEVELOPMENT OF MATHEMATICAL MODEL

OF A RADIAL SOLID METAL CORD TYRE WITH HYPERELASTIC

PROPERTIES OF RUBBER TAKEN INTO ACCOUNT

O.A. Odintsov

Bauman Moscow State Technical University, Moscow, Russian Federation e-mail: odintsov@bmstu.ru

A mathematical model of the automobile radial solid metal cord tyre (a model of geometrically nonlinear anisotropic three-layer shell) is presented which takes into account the physical nonlinearity of a matrix of cord-rubber composite layers. The determining relationships based on the two-term hyperelastic Mooney-Rivlin potential that was added with a penalty term for taking into account the low rubber compressibility were used for the rubber matrix simulation. Results of solving the test nonlinear problems of elasticity theory as well as the problems of quasi-static loading of a fast fixed radial tyre under conditions of the contact with a smooth flat supporting surface are given for two types of tyres (the motor-car tyre with a textile frame and the truck solid metal cord tyre).

Keywords: finite element method, hyperelasticity, rubber, tyre.

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

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

ISSN 0236-3941. Вестник МГТУ им. Н.Э. Баумана. Сер. “Машиностроение” 2013. № 3 11

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

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

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

Рис. 1. Основные элементы конструкции радиальной шины

12 ISSN 0236-3941. Вестник МГТУ им. Н.Э. Баумана. Сер. “Машиностроение”. 2013. №3

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

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

Определяющие соотношения для резины. Для гиперупругих материалов определяющие соотношения, как правило, формулируются функциями упругого потенциала. В случае изотропии упругих свойств материала, что в достаточной степени истинно для резины, упругий потенциал выражается через инварианты деформированного состояния, например через инварианты правого тензора Коши-Грина: C = F1F. Выражения для инвариантов тензора C имеют следующий вид:

3

Ic = tr C = ^ cn; i=1 33 (1)

IIc = tr(C : C) = ££ cj; i=1 j=1 (2)

IIIc = det C = (det F )2 = J2. (3)

Поскольку третий инвариант правого тензора Коши-Грина равен квадрату определителя J тензора градиента места F, эти два параметра взаимозаменяемые.

Одной из наиболее простых моделей, описывающих гиперупругое поведение резины, является модель Муни-Ривлина [1, гл. 6.6.2], имеющая следующий вид в двухконстантной формулировке:

W = W (Ic ,IIc ,J ) =

= ^io (Ic — 3) + ^^oi (i2 — IIc — 6) + ^к (j — 1)2 , (4)

где ^10, ^01 — постоянные материала; к — параметр штрафа, имеющий смысл модуля объемной деформации для гомогенизированного

упругого потенциала W = W\J-3Ic, J-3IIc, J). В случае малых

деформаций G = 2 (^10 + ^01), где G — модуль упругости 2-го рода.

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

ISSN 0236-3941. Вестник МГТУ им. Н.Э. Баумана. Сер. “Машиностроение” 2013. № 3 13

выражений для тензоров напряжений и деформаций, путем выполнения интегрирования по начальному объeму V:

где A — работа постоянных внешних сил на достигнутых перемещениях.

Тогда принцип возможных перемещений будет описываться вариационным уравнением вида

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

где n — размерность задачи.

Таким образом, все характеристики деформированного состояния также могут быть однозначно выражены через вектор узловых степеней свободы {q}; следовательно, функциями вектора {q} являются как инварианты тензора деформаций, так и сам упругий потенциал:

ic = Ic({q}); lie = lie({q}); J = J({q}); W = W({q}). (8)

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

Следовательно, 8W = (V{q}W){dq}.

lxn nxl

Это позволяет перейти от вариационного уравнения (6) к системе нелинейных алгебраических уравнений вида

где {F} — вектор узловой нагрузки; {R} — узловая невязка сил.

(5)

(6)

{и} = [N ]{q}

(7)

3x1 3xnnx1

(9)

Для краткости обозначим

(V(b}{a})

= da

j = дь,.

j

(10)

(11)

14 ISSN 0236-3941. Вестник МГТУ им. Н.Э. Баумана. Сер. “Машиностроение”. 2013. №3

Наиболее эффективным способом решения уравнения (11) является использование итерационной схемы, основанной на методе Ньютона-Рафсона или его модификациях. Для построения системы разрешающих уравнений на k-м шаге решения необходимо выполнить линеаризацию невязки {R} вблизи текущего приближения {qk}, что приводит к следующей системе линейных уравнений относительно искомого приращения узловых степеней свободы {Aqk}:

{Rk} + V{q}{R}{Aqk } = 0 (12)

nx 1 nxn nx 1

или

[Кт ]{A qk } = -{Rk},

nx n nx 1 nx 1

(13)

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

[К ] = V,S}{R} =

nx n n n

V

V {q}{q}WdV.

(14)

В выражении (14) оператор V{q}{q} W обозначает вычисление матрицы Гессе функции W({q}), а внешние силы принимаются не зависящими от перемещений.

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

Эту процедуру разобьем на три этапа. На первом этапе сформируем из компонент тензора градиента места F в глобальной декартовой системе координат арифметический вектор {0} размера 9 х 1:

m

9x1

ди1 ди1 ди1 ди2

1 +-------- ---- ----

дХ^ дХ^ дХ3 ’ дХ^

1 +

ди2 дХ2,

ди3 i + ди3 ..,дХ2, + дХз

т

,

(15)

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

{0} = {I} + [dN]{q}, (16)

9x 1 9x n 9x 1

где {I} — вектор постоянных коэффициентов, принимающих значения 0/1; [dN] — матрица производных функций формы, получаемая с использованием изопараметрической техники [2, гл. 9.2] и не зависящая от {q}, тогда

V{q}{0}

9x n

т

дqj

[dN] .

(17)

ISSN 0236-3941. Вестник МГТУ им. Н.Э. Баумана. Сер. “Машиностроение” 2013. № 3 15

Считая, что W = W({0({q})}), выразим первые и вторые производные упругого потенциала W в виде функции вектора {0}:

V(q}W = VWWV{q}{0} = VmW [dN]; (18)

1xn 1x9 9xn

V{,mW = (V{S}{#})T (V{,}{,}W )(V{,}{«}) =

nx n nx 9 9x 9 9x n

= [dN]T (V{o}{o}W) [dN]. (19)

Следующим шагом раскроем представление упругого потенциала W как

W({0}) = W({I({0})}), {I} = [1с({0}),На({0}),J({0})],т (20)

3x1

и выполним дифференцирование W как сложной функции:

VWW = V{J }W VW{I}; (21)

1x9 1x3 3x9

V,»H„W = (V,„,{I })T(V„№ }W )(V,»,{I}) + VmW VWW{I},

9x9 9x3 3x3 3x9 1x3 3x9x9

(22)

что с учетом (4) приводит к следующим выражениям:

V{<9}W = (^10 + ^01IC) V{d}IC — 2^01V{<9}IIC + к (J — 1) V{e}J; (23)

V{<9}{<9}W = V{d}Ic) T (V{0}1^ + (^10 + ^01IC) V{в}{в}1с —

9x1 1x9 9x9

— 2 ^01V{0}{0}IIC + K( V{d}J) T (V{<9} J) + K (J — 1) V{d}{e}J. (24)

2 9x9 9x1 1x9 9x9

Заключительным этапом построения матрицы тангенциальных жесткостей является получение значений градиента и матриц Гессе для инвариантов тензора деформаций как функций вектора производных {0}.

Для этого отметим, что из определения вектора {0} и тензора деформаций C непосредственно следует, что

Ic = tr C = tr ([F]T [F]) = [F] : [F] = {0}T {0} ; (25)

V{o}Io = 2{0}T; (26)

V{e}{e}IC = 2[I], (27)

9x 9

где [I] — единичная матрица.

Чтобы получить аналогичные выражения для IIC, выполним преобразование компонент тензора C в вектор-столбец подобно тому, как

16 ISSN 0236-3941. Вестник МГТУ им. Н.Э. Баумана. Сер. “Машиностроение”. 2013. №3

это было сделано для тензора градиента места в (15):

Тогда

{c}

9x1

( $2 + $2 + $2 ч

$1$4 + «2«5 + $3$6

$1$7 + $2$8 + $3$9 $1$4 + $2$5 + $3$6

< $4 + $2 + $2 >.

$4$7 + $5$8 + $6$9 $1$7 + $2$8 + $3$9 $4$7 + $5$8 + $6$9 $2 + $8 + $2

По = {c}T {c}

(28)

(29)

и с учeтом этого можно получить следующие выражения:

Vm lie = V{c) lie Vm{c}; (30)

1x9 1x9 9x 9

V{c) lie = 2{c} T. . (31)

Vm lie = 2{c}T Vm{c} (32)

где матрица V(^){c} имеет вид

"2$1 2$2 2$з 0 0 0 0 0 0 '

$4 $5 $6 $1 $2 $3 0 0 0

$7 $8 $9 0 0 0 $1 $2 $3

$4 $5 $6 $1 $2 $3 0 0 0

Vrn{c} = 0 0 0 2$4 2$5 2$6 0 0 0 . (33)

9x 9 0 0 0 $7 $8 $9 $4 $5 $6

$7 $8 $9 0 0 0 $1 $2 $3

0 0 0 $7 $8 $9 $4 $5 $6

0 0 0 0 0 0 2$7 2$8 2$9

Вторые производные инварианта 11с, учитывая (32), можно записать как

д211с =2sr- (dCkdCk + d2ck \ (34)

д$гд$3 = \d$j д$г + Ckд$гd$J ( )

или в матричной форме как:

Ис = 2(V,„{c})T(VM{c>) +4 [A], (35)

9x 9 9x 9

где симметричная матрица [A] составлена из компонент вектора {c}

ISSN 0236-3941. Вестник МГТУ им. Н.Э. Баумана. Сер. “Машиностроение” 2013. № 3 17

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

С1 0 0 С2 0 0 С3 0 0

0 С1 0 0 С2 0 0 С3 0

0 0 С1 0 0 С2 0 0 С3

С4 0 0 с5 0 0 Сб 0 0

A = 0 С4 0 0 с5 0 0 С6 0 . (36)

9x9 0 0 С4 0 0 с5 0 0 С6

С7 0 0 С8 0 0 С9 0 0

0 С7 0 0 С8 0 0 С9 0

0 0 С7 0 0 С8 0 0 С9

Аналогичные выражения для определителя J могут быть получены непосредственным дифференцированием:

J

V{^} J

016569 — 0\0б0^ — 646269 + 646368 + #7#2#6 6569 — в 6 6 8

вгвб — 64 69

в4 в8 — 6765 6368 — 62 69

< 6169 — 6763 > ;

6762 — 6168

6266 — 63 65 6463 — 616б к 6165 — 64 62 ,

676365; (37)

(38)

0 0 0 0 69 — 68 0 — 6б 65

0 0 0 — 69 0 67 66 0 — 64

0 0 0 68 —67 0 — 65 64 0

0 —69 68 0 0 0 0 63 —62

VmmJ = 69 0 — 67 0 0 0 —63 0 61

— 68 67 0 0 0 0 62 —61 0

0 66 — 65 0 —63 62 0 0 0

—6б 0 64 63 0 —61 0 0 0

65 —64 0 — 62 61 0 0 0 0

Итоговые выражения для градиента и матрицы Гессе упругого потенциала как функции узловых степеней свободы {q} получаются путем подстановки выражений (25)-(27), (29), (32), (35) и (37)-(39) в выражения (23) и (24) соответственно, с последующим преобразованием результата с использованием выражений (18) и (19). Интегрирование по объему выполняется обычным для МКЭ способом с использованием сетки конечных элементов и квадратурных формул Гаусса-

18 ISSN 0236-3941. Вестник МГТУ им. Н.Э. Баумана. Сер. “Машиностроение”. 2013. №3

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

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

тензора Коши-Грина: E = - (C — I). В качестве тензора напряжений

2

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

о =

2 dW. 4 dW

Ь ^77-Ь2 + pI,

(40)

JdIc J dIIc

где b = FFT — тензор Фингера, p = k(J — 1) — гидростатическое давление.

Эффект объемного заклинивания. При непосредственном использовании дискретизированной формы уравнения (6) возникает хорошо известный эффект “объемного заклинивания” [3, гл. 8], проявляющийся в значительном ужесточении конструкции в условиях стесненных деформаций и являющийся следствием несогласованности порядка аппроксимации перемещений и гидростатического давле-

dW (J 1)

ния p = J = k(J — 1).

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

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

ISSN 0236-3941. Вестник МГТУ им. Н.Э. Баумана. Сер. “Машиностроение” 2013. № 3 19

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

Разобьeм гомогенизированный упругий потенциал W на два слагаемых, W0 и Wp, первое из которых связано с энергией формоизменения, а второе формирует энергию объeмной деформации:

W = Wo (Ic ,IIc ,J) + WP(J); (41)

Wo(Ic ,IIc ,J ) =

= ^10 J 3 Ic — ^^01 ^J 3 (IC — IIC — ^ ; (42)

1 9

Wp(J) = 2К (J - 1)9. (43)

Выполним раздельное построение матрицы [KT]:

no no no

[k ] = £££ giigi2giз det V{£}{X}V{q}{q} ^^o+

nn

i1 = 1 %2 = 1 i3 = 1

Пр Пр Пр

+ EEE gilgi2gi3 det V{?}{X}^^{q}{q}^^p, (44)

il = 1 i 2 = 1 i3 = 1

где n0, np — порядок схемы интегрирования энергии формоизменения и энергии объeмной деформации соответственно; gi — весовые коэффициенты квадратурной формулы Гаусса-Лежандра соответствующего порядка, определитель det V{g}{X} — якобиан изопараметрического преобразования координат.

Аналогичным образом выполняется построение вектора-столбца невязки узловых сил {R}.

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

трилинейных функций формы Щ£,г},С) = -(1 + &£)(1 + цм)(1 + CiC)

8

необходимо выбирать np = 1, тогда как параметр no для недопущения возникновения механизмов на элементе должен удовлетворять условию n0 ^ 2. В [3, гл. 8] указано, что подобная схема построения элемента в ряде случаев может давать ложные осцилляции гидростатического давления, поэтому для практического применения описанного подхода необходимо выполнить тестовые расчеты.

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

Первая задача заключалась в расчете плоского деформированного состояния клиновидного уплотнителя (рис. 2), нагруженного постоян-

20 ISSN 0236-3941. Вестник МГТУ им. Н.Э. Баумана. Сер. “Машиностроение”. 2013. №3

Рис. 2. Исходная и деформированная конфигурации уплотнителя в Ansys

ным давлением по верхней грани и жестко заделанного снизу. Размеры сечения: высота 3 мм, ширина снизу 1 мм, ширина сверху 2 мм. Характеристики материала: д10 = 1,071 МПа, ^01 = 0, к = 6000 МПа (неогуков материал). Давление составляло: 1,5 МПа. При решении в Ansys использовался элемент PLANE182 с билинейными функциями формы.

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

Аналогичная картина наблюдается и для графика распределения гидростатического давления в том же сечении (рис. 4). В то время как схема сокращенного интегрирования дает решение, полностью совпадающее с решением в комплексе Ansys, применение полного интегрирования существенно искажает решение.

Рис. 3. Распределение прогиба в центральном сечении уплотнителя

ISSN 0236-3941. Вестник МГТУ им. Н.Э. Баумана. Сер. “Машиностроение” 2013. №3 21

Рис. 4. Распределение гидростатического давления в центральном сечении уплотнителя

Также была решена объемная задача растяжения стержня прямоугольного сечения из гиперупругого материала (рис. 5), торцы которого были жeстко заделаны и квазистатически нагружены силами, изменяющимися от нуля до заданного значения. Размеры стержня: длина 100 мм, высота и ширина сечения 15 мм. Характеристики материала: Мю = 0,50783 МПа, д01 = 0,54114 МПа, к = 6000 МПа; сила составляла 1500 Н. При решении в Ansys использовался элемент SOLID185 с трилинейными функциями формы. Число элементов в каждом из направлений: 5 х 5 х 34. Дополнительно было выполнено решение в комплексе Abaqus с использованием элементов C3D8.

Сравнение нагрузочных характеристик для стержня (рис. 6) показывает, что так же, как и для плоской задачи, в объeмном случае схема выборочного интегрирования дает решение, совпадающее с решением в Ansys и Abaqus. Также было обнаружено, что при выборе для решения этой задачи в Abaqus конечного элемента C3D8R, использующего сокращeнное интегрирование, решение начинает расходиться

Рис. 5. Исходная и деформированная конфигурации стержня в Ansys

22 ISSN 0236-3941. Вестник МГТУ им. Н.Э. Баумана. Сер. “Машиностроение”. 2013. №3

100 150 200 250 300 350 400 450 500 550 Удлинение, %

Рис. 6. Нагрузочная характеристика стержня

при нагрузке, составляющей примерно 70% заданной, при этом для меньших значений нагрузки такое решение соответствует решению, получаемому с использованием элементов C3D8.

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

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

Вычислительные эксперименты показали, что само по себе использование квадратичных функций формы значительно снижает влияние эффекта объемного заклинивания по сравнению с использованием линейных функций (с примерно 40 % до 5 %), однако это влияние всe равно нельзя игнорировать, поэтому в 18-узловом элементе также была использована схема выборочного интегрирования. Устойчивая

ISSN 0236-3941. Вестник МГТУ им. Н.Э. Баумана. Сер. “Машиностроение” 2013. №3 23

Рис. 7. Схема 18-узлового конечного элемента трехслойной оболочки

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

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

Для оценки корректности результатов расчeта, их сравнивают с результатами, полученными с использованием модели трсхслойной оболочки с линейной прослойкой [7]. Использование физического эксперимента для подобного сопоставления затруднено, поскольку на практике силы в корде экспериментально не измеряются, а влияние выбора модели прослойки на интегральные характеристики шины, измеряемые в эксперименте, незначительно. Сопоставление проводилось на двух типах шин: легковой радиальной шине типоразмера 175/70R13, имеющей однослойный текстильный каркас и два слоя металлокорда в брекере, и грузовой радиальной ЦМК-шине типоразмера 295/80R22,5, у которой в каркасе используется один металлокордный слой, а бре-керный пояс состоит из четырех перекрeстно армированных металлокордных слосв.

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

24 ISSN 0236-3941. Вестник МГТУ им. Н.Э. Баумана. Сер. “Машиностроение”. 2013. №3

Рис. 8. Результаты расчета для шины 175/70R13

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

При расчeте ЦМК-шины по модели трeхслойной оболочки с линейно-упругой прослойкой наблюдается как ухудшение скорости сходимости итерационного решения, так и возникновение в корде областей резкого изменения сил без наличия каких-либо объективных причин для этого. На рис. 9 показано распределение сил в корде двух наиболее нагруженных слоeв брекера грузовой ЦМК-шины. На графике видно, что в третьем слое брекера (B3, нумерация с внешней стороны) вблизи центра пятна контакта наблюдаются всплески сил, хотя и геометрия шины, и распределение нагрузки в этой области являются гладкими, а геометрические и механические свойства внутренних слоeв шины изменяются плавно.

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

Рис. 9. Результаты расчета для шины 295/80R22,5 с линейно-упругой прослойкой

ISSN 0236-3941. Вестник МГТУ им. Н.Э. Баумана. Сер. “Машиностроение” 2013. №3 25

Рис. 10. Результаты расчета для шины 295/80R22,5 с гиперупругой прослойкой

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

Также было получено, что несмотря на увеличение времени построения матрицы тангенциальных жесткостей на отдельной итерации (с 1,15 до 1,61 с) за счeт применения более сложной модели прослойки, общее время вычислений в итоге значительно сокращается благодаря лучшей сходимости решения на отдельных шагах нагружения (253 с против 342 с на расчeт).

Заключение. В работе представлена модель трeхслойной геометрически нелинейной оболочки, учитывающая физическую нелинейность резиновой прослойки между каркасом и брекером посредством использования двухконстантного гиперупругого потенциала Муни-Ривлина. Был построен комбинированный объемно-оболочечный конечный элемент, имеющий функции формы пониженного порядка по направлению толщины прослойки и использующий метод выборочного интегрирования для устранения эффекта объемного заклинивания.

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

ЛИТЕРАТУРА

1. Bathe К.-J. Finite element procedures. Prentice Hall, 1996.

26 ISSN 0236-3941. Вестник МГТУ им. Н.Э. Баумана. Сер. “Машиностроение”. 2013. №3

2. Zienkiewicz O.C., Taylor R.L. The finite element method. Oxford: Butterworth-Heinemann, 2000. Vol. 1.

3. Belytschko Т, Liu W.K., Moran B. Nonlinear finite elements for continua and structures. New Jersey, John Wiley Sons Inc. 2000.

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

5. Белкин А.Е., Одинцов О.А. Численное решение геометрически нелинейной задачи контакта автомобильной шины с твердой опорной поверхностью // Вестник МГТУ им. Н.Э. Баумана. Сер. Машиностроение. 2007. № 1 (66). С. 20-35.

6. Белкин А.Е., Одинцов О.А. Численное решение контактной задачи для автомобильной шипы при стационарном качении с боковым уводом // Изв. вузов. Машиностроение. 2008. № 6. С. 27-42.

7. Одинцов О.А. Численный анализ контакта автомобильной шины в конечноэлементном комплексе Cascad FEA // Проблемы шин и резинокордных композитов: Труды XXI Междунар. конф. М., 2010. Т. 2. С. 101-107.

REFERENCES

1. Bathe K.J. Finite Element Procedures. New Jersey, Prentice Hall, 1996.

2. Zienkiewicz O.C., Taylor R.L. The finite element method. Oxford, Butterworth-Heinemann, 2000, vol. 1.

3. Belytschko T, Liu W.K., Moran B. Nonlinear finite elements for continua and structures. New Jersey, John Wiley & Sons, 2000. 641 p.

4. Belkin A.E. Razrabotka sistemy modeley i metodov rascheta napryazhenno-deformirovannogo i teplovogo sostoyaniy avtomobil’nykh radial’nykh shin. Diss. dokt. tekhn. nauk [The development of models and methods for calculating stress-strain and thermal states of automotive radial tires. Dr. tech. sci. diss.]. Moscow, MGTU im. N.E. Baumana Publ., 1998. 284 p.

5. Belkin A.E., Odintsov O.A. A numerical solution to the geometrically nonlinear problem for a tire contact with a solid support surface. Vestn. Mosk. Gos. Tekh. Univ. im. N.E. Baumana, Mashinostr. [Herald of the Bauman Moscow State Tech. Univ., Mech. Eng.], 2007, vol. 66, no. 1, pp. 20-35 (in Russ.).

6. Belkin A.E., Odintsov O.A. A numerical solution to the contact problem for a car tire under steady-state rolling with side slipping. Izv. Vyssh. Uchebn. Zaved., Mashinostr. [Bull. Inst. Higher Educ., Mach. Constr.], 2008, no. 6, pp. 27-42 (in Russ.).

7. Odintsov O.A. Numerical analysis of the car tire contact in finite element package “Cascad FEA”. Tr. 21 Mezhdunar. Konf. “Problemy shin i rezinokordnykh kompozitov” [Proc. 21st Int. Conf. “Problems of tires and rubber-cord composites”]. Moscow, 2010, vol. 2, pp. 101-107 (in Russ).

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

Олег Александрович Одинцов — доцент кафедры “Динамика и прочность машин” МГТУ им. Н.Э. Баумана. Автор 12 научных работ в области контактных задач механики деформируемого твердого тела, вычислительных методов нелинейной теории упругости.

МГТУ им. Н.Э. Баумана, Российская Федерация, 105005, Москва, 2-я Бауманская ул., д. 5.

O.A. Odintsov — assoc. professor of “Dynamics and Strength of Machines” department of the Bauman Moscow State Technical University. Author of 12 publications in the field of contact problems of mechanics of deformable solids, computational methods of nonlinear theory of elasticity.

Bauman Moscow State Technical University, Vtoraya Baumanskaya ul., 5, Moscow, 105005 Russian Federation.

ISSN 0236-3941. Вестник МГТУ им. Н.Э. Баумана. Сер. “Машиностроение” 2013. № 3 27