К вопросу расчетной оценки несущей способности кузовных конструкций Текст научной статьи по специальности «Физика»

УДК 629.113

Е.А. Наумов

К ВОПРОСУ РАСЧЕТНОЙ ОЦЕНКИ НЕСУЩЕЙ СПОСОБНОСТИ КУЗОВНЫХ КОНСТРУКЦИЙ

Нижегородский государственный технический университет им. Р.Е. Алексеева

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

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

Во многих прикладных отраслях строительной механики важна задача, связанная с оценкой несущей способности инженерных конструкций. Кузовостроение не является исключением, где несущая способность является важной частью комплексного понятия пассивной безопасности, которое, в свою очередь, регулируется нормативно-правовыми документами. Актуальность проблемы повышения пассивной безопасности автомобилей и автобусов [2, 3] требует постоянного совершенствования методов ее оценки на всех стадиях проектирования. При этом на начальных этапах проектирования кузовных конструкций, при выборе их силовых схем, если позволяет конструкция, целесообразно использование балочных моделей. На рис. 1 показано нагружение автомобиля при испытаниях по методике EuroNCAP на боковое столкновение с единичным препятствием (pole-тест) и балочная схематизация салона кузова легкового автомобиля при рассмотрении нагружения центральной стойки одной из составляющих боковой аварийной нагрузки.

Рис. 1. Схема испытаний по EuroNCAP Pole Protocol и расчетная модель салона кузова легкового автомобиля

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

© Наумов Е.А., 2014.

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

На рис. 2 показан фрагмент конструкции, состоящий из двух элементов (I и к) с общим узлом N. Элементы I и к в узлах Р и М соединяются с другими элементами.

Рис. 2. Фрагмент КЭМ конструкции

В разработанном алгоритме за основу математической модели балочного элемента взята формулировка элемента Белычко (Belytschko). В ней используется метод вращения координат (co-rotational technique) для обработки больших углов поворота. На рис. 3 изображен изопараметрический балочный элемент [6]; тройка единичных векторов {ei} образует элементную систему координат (element coordinate system). Орт e1 всегда направлен вдоль оси элемента от начального (I) узла к конечному (J).

Рис. 3. Балочный конечный элемент (параметры 8 - геометрические характеристики сечения)

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

т = т = «А0; Л 8ЛЛ. (1)

где т - масса; Jik - компоненты тензора инерции; р - плотность материала; А - площадь по-

перечного сечения; 10 - длина элемента в недеформированном состоянии; - единичный тензор; Jk - момент инерции сечения относительно локальной оси с ортом ек.

Здесь и далее при описании векторных и тензорных величин используется тензорная нотация, когда в соответствующих выражениях индекс, встречающийся раз, принимает значения 1, 2, 3, а по дважды повторяющемуся индексу в одночлене подразумевается суммирование от 1 до 3. Символ «Л» над величинами означает, что они определены в локальной системе координат элемента |е;}.

В ЭСК тензор инерции диагоналей, т.е. оси ЭСК являются его главными осями.

Суммарная масса узла N на рис. 2 определяется как сумма узловых масс элементов I и к. Для нахождения суммарного узлового тензора инерции необходимо преобразовать компоненты тензоров элементов в одну и ту же систему координат. В данном случае используется глобальная систему координат (ГСК). Для тензора 2-го ранга в декартовой прямоугольной системе координат преобразование компонент при смене базиса [1] имеет вид

•1¡к = еИект^1т, (2)

где к - компоненты тензора в новом базисе (в данном случае в ЭСК); Jlm - компоненты

тензора в старом базисе (в данном случае в ГСК); ек - косинус угла между ьй осью ЭСК и к-й осью ГСК (или координата к вектора е; в ГСК). Обратное преобразование соответственно

¿¡к = еИетк^ 1т. (3)

Преобразование (3) эквивалентно в матричном виде можно записать как

J = ^ 1 , (4)

где ^ - матрица преобразования, столбцы которой - векторы ei.

Итак, после преобразования тензоров инерции элементов I и к суммарный тензор инерции для узла N будет вычислен как их сумма.

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

а(Ю = м-1^) -$) -гОО); у(п+1/2) = у(п-1/2) + а(й)Л^+1/2);

+ л/и+1) (5)

и(п+1) = и(п) + у(«+1/2)л^(«+1). ЛГ(И+1/2) = Л + Л

; 2 ' где (5) М - матрица масс; - вектор внешних приложенных сил; {-(¡1 - вектор внутренних сил в узле (вычисляется суммированием векторов внутренних усилий по всем элементам, соединяющимся в данном узле); {- вектор сил вязкого сопротивления; Дt - шаг интегрирования. Матрица масс узла определяется суммарной узловой массой (компоненты, соответствующие поступательным степеням свободы) и суммарным тензором инерции (компоненты, соответствующие вращательным степеням свободы).

Очевидно, что уравнения (5) наиболее просто (в разделенном виде) решаются в случае с диагональной матрицей масс. Масса - это скаляр, который не меняется при смене базиса. Поэтому для поступательных степеней свободы никаких преобразований выполнять не нужно. Иначе обстоит дело с угловыми степенями свободы. Выше указывалось, что тензор инерции узлов элемента в ЭСК задан своими главными значениями. Однако суммарный узловой тензор инерции в общем виде не является диагональным, поэтому в ГСК решить уравнения, связанные с определением ускорений, в разделенном виде не получится. Указанное затруднение решается через проблему собственных значений для суммарного узлового тензора инерции, когда определяются направления главных осей инерции, в которых центробежные компоненты тензора равны нулю, а на главной диагонали стоят его главные значения. На рис. 2 направления главных осей инерции в узлах показаны тройками ортов {Ь;}, определяющих узловые системы координат (УСК). УСК вращается вместе с узлом.

Необходимо отметить, что при вычислении ускорений вращательных степеней свободы по формуле (5) вектор суммарного узлового момента (выражение в скобках) должен быть преобразован в УСК. Компоненты вектора Г при изменении системы координат преобразуются как [1]

Ц = ъ.кЦк, Ц = ъкгЦк, (6)

где Ц - компоненты вектора в УСК; Ъ± - косинус угла между 1-й осью УСК и к-й осью ГСК. Или эквивалентным образом в матричном виде

Ж = X Ё; Ё = ХТ Ж, (7)

где X - матрица преобразования, столбцы которой векторы Ьг-.

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

£Т = \/х1 > тхЗ > т у1 > ту1> тг1> } (8)

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

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

£0(п+1) = £(п) + (п) (9)

После этого осуществляется проверка на достижение предельного состояния в каждом узле элемента. Если предельное состояние не было достигнуто, то полагается

£(п+1) = £0(п+1) и осуществляется переход к следующему шагу. В противном случае осуществляется специальное корректирование вектора усилий, вычисленных по формуле (9). Подробнее об этом далее в статье.

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

% = I -1° . (10)

Изгибные деформации определяются разложение векторного произведения текущего орта е1 на его исходное направление, поворачивающееся вместе с узлом, в текущем базисе |е;}:

е° 4°. }= ^} пп

е1 X е° = 0уе2 + 02еЪ = -е?13е2 + е1 2е3 •

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

вхЛ =1 е1 (е2! х еЪ1 + еЪ1 х е° )=1 (е?21е°21 - е?21ей1 + еУз1е°31 - е?31е?31). (12)

В итоге, вектор деформаций

¿Т = &хи> 0у!, 6и, вг1, }. (13)

В формулах (10)-(12) в качестве исходного взято состояние в начальный момент времени (¿0), поэтому вычисляемые деформации являются суммарными. Чтобы получить приращения деформаций, нужно в качестве исходного взять состояние на предыдущем шаге, а в формулах заменить 0 (нуль) на (п-1), полагая, что текущая конфигурация системы рассматривается в момент времени /п) (в выражениях (10)-(13) индекс п, соответствующий текущему состоянию, опущен. При этом суммарные деформации:

i(«+i) _ А(")

= d(n) + Ad'

(n)

(14)

Вернемся к вопросу о вычислении внутренних усилий. Итак, после вычисления приращения деформаций (см. формулы (10)—(12)), рассчитываются векторы приращений и суммарных внутренних усилий по формуле (9). Отметим, что остальные внутренние силовые факторы определяются из условия равновесия:

jxl = -f—; mxl = -mxJ; fzl = "

myl + myJ

; f zJ = -f zJ; f yJ =--z^+——; f yl = -fyJ • (15)

I I

Далее предлагается использовать следующее условие предельного состояния в узле:

h(f ) = 1

2

kv

J V

J x Fx0

f

+

—

y

f Г

M

y0

+

—,

Л

2

v mz0j

+ k„

f Л2

Jy

v Fy0

+

Fz0

2

f

+

—

\

2

v Mx0 j

-1

< 0,

(16)

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

Fx0 = Aar>'

My,z0 = Wrrn y, z°T;

Fy , z0 = К zAtt;

y,z

Mx0 = Wrn x Тт ,

(17)

где A - площадь поперечного сечения; Wпл - пластические моменты сопротивления относительно соответствующих осей; от - предел текучести при растяжении; тт - предел текучести на сдвиг; ky,z - коэффициенты сдвига в направлении локальных осей у и г соответственно.

После вычисления усилий по формуле (9) для каждого узла балочного элемента проверяется условие (16). Если оно нарушается, то задача заключается в корректировке усилий таким образом, чтобы неравенство к > 0 обратить в равенство. Для упрощения выражений введем следующие обозначения, опуская индекс п для номера шага по времени:

{¿/>туЬ т21> туЗ' mzJ }= {Х }

Ко, Мх0, Му0, Ы20 }= X, X2, X3, X4 }■ (18)

Х3 + Х5 (х6 = Ру о1 ) = - ^ X = Fzоl )

f zl =-"

I V 6 у ^ I

Выражение (16) для узлов I и J с учетом введенных обозначений принимает вид

hi(x )=2 h2 (x) = 1

\2

\2

(

\2

Л

V

x1 + x2 + x3 + x4 + x4 + x6 + x3 + x5

V Xi j V X 2 J V X 3 J V X 4 J V X5 J V X 6 J

2

x1 + x2 + x5 + x6 + x4 + x6 + x3 + x5

V Xi j V X 2 J V X 3 J V X 4 J V X 5 J V X 6 J

-1

-1

< 0;

< 0.

(19)

При нарушении условий (19) возможны два случая.

Случай 1 - условие (19) нарушено в одном узле, а в другом выполнено.

Для примера рассмотрим ситуацию, когда нарушено условие к1, а к2 выполнено. Рассуждения в обратной ситуации проводятся аналогично. Итак, при к1(х0) = а > 0 необходимо произвести корректировку усилий в 1-м узле (I). При этом усилия (изгибные моменты) от второго узла, входящие в выражение, не корректируются. Возврат усилий на предельную поверхность осуществляется их умножением на корректирующий скаляр в, определяемый как решение квадратного уравнения с коэффициентами а, Ь, с:

2

2

2

2

2

2

2

2

а = Х

г=1

( „О Л

( „0 V

-г + —4 + ; Ъ = 2

V V V X 0 у VX 6у

—л X,

4 Х6

V ^5

2

+ -

Х3 Х5 ^

X

2

О 2 О 2

С =

—6 + —0

V X0 у V X у

-1; Б = Ъ2 - 4ас; в =

6 у

-Ъ+УБ

2а

(20)

х =

вО вхО, вХз0, в*0 —0 хО }

Случай 2 - условие (19) нарушено в обоих узлах.

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

Ь, (х° )= Ь2 (х0 )= 1

( о V

XI

+

( 0^2 (у 0 ^ Х3(0)

V Х2 у

+

X,

(„0 V

+

х,

4(6)

X.

-1

= а;

(21)

вх=Рх°.

л/1 + 2а

Второй, более сложный вариант, будет, когда предельные условия нарушены в обоих узлах и при этом не равны между собой (М(х°) Ф Л2(х°)). В этом случае предлагается вектор внутренних усилий в предельном состоянии (х) искать как решение задачи нелинейного программирования о поиске минимального расстояния в относительных координатах, которым алгоритм пластичности будет переводить точку из начального состояния (х°) в искомое (формулы (22) ниже). Условия И(х) = ° соответствуют нахождению точки искомого состояния на линии пересечения предельных поверхностей.

Точка х должна удовлетворять условиям И1(х) = И2(х) = ° с учетом согласованности знаков скорректированных ВСФ и деформаций, которые их вызвали (для этого и нужны ограничения в виде неравенств).

1 6

Г (х ) =1X

0_ >

2 /=1

Н

^ шт, Нг =

, х0 * 0 _ , г = 1,6;

1, хО = 0

Ь,« = 2 И2 (х) = 1

—1 + —2 + —3 + —4 + —4 + —6 + —3 + —0

V V X 2 у V X 3 у V X 4 у V X0 у V X6 у

V V г „ л

2

2

—1 + —2 + —0 + —6 + —4 + —6 + —3 + —0

V у V X 2 у V X 3 у V X 4 у V X0 у V X6 у

-1

-1

(0х0).

/ ч — - шах10,—■ I п ~ / ч

ёг (х) = -—^^ < 0; ёг (х) =

Н

шт (0,—)0)- —

= 0;

= 0;

(22)

Н

< 0, г = 1,6.

Задача (22) может рассматриваться и как альтернатива нахождению вектора предельных усилий х с помощью масштабирующего множителя в по формулам (2°)-(21). При этом из условий-ограничений исключены неравенства и одно равенство (в рассмотренных примерах к2), а размерность задачи уменьшена до четырех.

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

2

2

2

2

0

х

<

2

2

2

2

2

2

2

2

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

Библиографический список

1. Борисенко, А.И. Векторный анализ и начала тензорного исчисления / А.И. Борисенко, И.Е. Тарапов. - 3-е изд. - М.: Высш. шк., 1966. - 252 с.

2. Наумов, Е.А. К расчетной оценке пассивной безопасности кузова легкового автомобиля при боковом столкновении / Е.А. Наумов, Л.Н. Орлов // Журнал автомобильных инженеров. 2009. № 5(58). С. 19-21.

3. Повышение безопасности кузовов автобусов на стадии проектирования / Л.Н. Орлов [и др.] // Проетирование, испытание, эксплуатация и маркетинг автотракторной техники: сб. научн. трудов к 60-летию каф. «Автомобили и тракторы». - Н. Новгород: НГТУ, 1997. - С. 230-238.

4. Тумасов, А.В. Оценка несущей способности каркаса кузова автобуса по результатам компьютерного моделирования / А.В. Тумасов, Л.Н. Орлов, П.С. Рогов // Труды Нижегородского государственного технического университета им. Р.Е. Алексеева. 2011. Т 78. № 3 (96). С. 150-156.

5. Tumasov, A. Analysis of Simulation Approaches for Calculative Estimation of Buses Passive Safety / A. Tumasov, A. Vashurin, P. Rogov, L. Orlov // Advanced Technologies for Commercial Vehicles - VTI2012E013: рroceedings of 2012 International Conference on Advanced Vehicle Technologies and Integration (VTI2012). 2012, Changchun, China. Р. 519-523.

6. Bathe, K.-J. Finite Element Procedures / Klaus-Jürgen Bathe // USA, Prentice-Hall, 1976. - 1050 p.

7. LS-DYNA THEORY MANUAL, March 2006, Compiled By John O. Hallquist.

Дата поступления в редакцию 22.04.2014

E.A. Naumov

THE ESTIMATED VALUE CAR AND BUS BODY SUPPORTING ABILITY

Nizhny Novgorod state technical university n.a. R.E. Alexeev

Purpose: Refined technique of estimation supporting ability of bodies during choosing the supporting scheme. Design/methodology/approach: Finite element method with explicit time integration scheme; co-rotational technique; limit state theory; tensor and vector algebra; nonlinear programming.

Findings: It is possible to apply the refined technique in perspective for more accurate estimation of supporting ability on the basis of solution simple (beam) models results that must shows higher convergence with experimental data. Research limitations/implications: The present paper provides a starting-point for further research in the field of structure supporting ability particularly automobile and bus bodies. In perspective It's expected to use plastic section modulus versus plastic strain relationship.

Originality/value: The main peculiarity of the study is modified limiting state criterion with 6th internal forces and technique to correct elastic prediction section forces increment based on conditional minimization (when plastic conditions both nodes exceed and unequal) that could provide more accurate values limit strength and energy.

Key words: limiting state criterion, body, supporting ability.