Конечный элемент в виде витка для расчета частот и форм собственных колебаний пружин c учетом сдвига и инерции вращения

Сорокин Ф.Д.; Су Чжоу

Расчет и конструирование машин

УДК 534.11:51-74 DOI 10.18698/0536-1044-2018-1-3-11

Конечный элемент в виде витка для расчета частот и форм собственных колебаний пружин c учетом сдвига и инерции вращения

Ф.Д. Сорокин, Су Чжоу

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

A Finite Element in the Form of a Single Coil for Calculating Natural Frequencies and Mode Shapes of Springs with Consideration of the Rotational Inertia and Shear

F.D. Sorokin, Su Zhou

BMSTU, 105005, Moscow, Russian Federation, 2nd Baumanskaya St., Bldg. 5, Block 1 Г(П] e-mail: sorokin_fd@mail.ru, szleon602@gmail.com

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

Ключевые слова: винтовые цилиндрические пружины, конечный элемент, виток, матрица жесткости, матрица масс, частоты и формы колебаний

To calculate natural frequencies and mode shapes of «long» cylindrical helical springs, a finite element in the form of a single coil is proposed, in which nodes are located on the axis of the spring. The stiffness matrix of the finite element is obtained by numerical integration of a system of differential equations for a spatial rod, therefore it can be regarded as numeri-

cally exact. The shape functions corresponding to the identity nodal displacements of the finite element are a by-product of the stiffness matrix calculation. The shape functions are also numerically exact as they are derived from a system of differential equations. The shape functions are used to build the mass matrix using the standard finite element method. When compiling the stiffness matrix, transverse shear and extension of the helical rod axis are taken into account, while the rotational inertia is considered when constructing the mass matrix. The main advantage of the developed finite element is that it is «soft» compared to a conventional beam finite element, whose tensile and bending stiffness values differ by orders of magnitude. In addition, the departure from partitioning the coil into rectilinear finite elements makes it possible to reduce the dimensionality of the problem by dozens of times. The comparison of the calculation and the experimental results has shown high accuracy of the proposed finite element.

Keywords: coil cylindrical springs, finite elements, stiffness matrix, mass matrix, frequencies and mode shapes

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

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

го витка (КЭВ) имеет определенные преимущества [2].

Цель работы — дальнейшее совершенствование КЭВ путем учета растяжения оси, поперечных сдвигов и инерции вращения, т. е. дополнительных эффектов, которые в работе [2] не учитывались.

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

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

Рис. 1. КЭВ цилиндрической пружины: Ра1-Раз и Рв1-Рвз — узловые силы в точках А и В; Ма1-Маз и Мв1-Мвз — моменты в точках А и В

никновению трудно контролируемой погрешности, особенно при вычислении старших форм колебаний [3-8].

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

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

Построение матрицы жесткости КЭВ. Для

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

dP п dM п п

— = 0; -+1X P = 0;

ds ds

^ = $X t + B • P; d*= A • M,

ds ds

(1)

где P, M — векторы внутренних сил и моментов в сечении; 5 — дуговая координата; t — орт касательной к оси стержня; u и Ф — векторы перемещений и поворотов сечения стержня; B — тензор податливости сечения от растяжения и сдвига; A — тензор податливости сечения от изгиба и кручения; X и • — знаки векторного и скалярного умножения.

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

Винтовую линию (ось пружины) задавали в обычном параметрическом виде:

r(s) =

( Ш^ *2(s) *s(s)

-R sin

s sin у s cos у

R cos

R

s cos у R

где у — угол подъема винтовой линии; R — радиус пружины (см. рис. 1).

Орт касательной определяли как производную от радиус-вектора оси:

dr t = —.

ds

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

A = — tt +—(E - tt); GJp EJX

„ 1 k ^ ч

B =-tt +-(E - tt),

EA GA

(2)

где G и E — модули упругости второго и первого рода; Jp и Jx — полярный и осевой моменты инерции сечения; E — единичный тензор; tt — диадное произведение касательного орта на себя (соответствует матричному умножению матрицы-столбца на матрицу-строку); A — площадь сечения; k — коэффициент учета неравномерности распределения касательных напряжений [15, 16] (для кругового сечения k = 10/9).

Площадь и моменты инерции круглого сечения определяли по формулам:

A=^dl т = т =^d! т = Kd4 A = 4 ; Jx =Jy = 64 ; Jp = 32 ,

где d — диаметр витка (см. рис. 1).

Элементы матрицы жесткости [K] КЭВ размерностью 12x12 находили последовательным решением двенадцати линейных краевых задач для системы дифференциальных уравнений (1). Численное решение каждой краевой задачи осуществлялось с помощью процедуры NDSolve программного обеспечения Wolfram Mathema-tica [17]. При этом, как и в работе [18], предвари-

тельно выполнялось преобразование векторных уравнений (1) к двенадцати уравнениям в проекциях на неподвижные декартовы оси хь x2, x3 процедурой Thread из того же программного обеспечения.

Вектор перемещений КЭВ размерностью 12x1 представляет собой объединение проекций узловых перемещений и поворотов на координатные оси:

Pa = —P(0); Ma = -M(0) - Re x P(0);

Y = (uA, «A, uB, «B)T,

(3)

где индексы «А» и «В» помечают соответствующие узлы (см. рис. 1), а направленные отрезки рассматриваются как матрицы-столбцы.

Для построения одного из столбцов матрицы жесткости в каждой из краевых задач одиннадцать узловых перемещений в выражении (3) задавались равными нулю, а узловое перемещение с номером i (г = 1, ..., 12) — единице. Например, для построения 4-го столбца матрицы [К]

Y = (0 0 0 1 0 0 0 0 0 0 0 0)Т. (4)

С помощью уравнений кинематики перемещения повороты из узлов КЭВ А и В переносились на края винтовой линии А0 и В0 (см. рис. 1), что позволило записать граничные условия для системы (1) следующим образом:

u(0) = uA +«A x eR; «(0) = «A; u(l) = ub +«в xeR; «(l) = «в,

(5)

где e — орт оси Хз; I — полная длина винтовой линии (оси) КЭВ, I = 2лЕ/со8 у.

Например, для вектора узловых перемещений, заданного соотношением (4), краевые условия для системы (1) имеют вид

(

u(0) =

0

—1-R

0

\

u(l) =

( 0 ^

0 0

(1 >

«(0) =

v 0 у 0

(6)

«(l) =

v 0 у

(7)

Решение линейной краевой задачи для системы (1) с граничными условиями (5) позволяет найти все компоненты напряженно-деформированного состояния винтового стержня, в том числе и силы, действующие на краях витка. Перенос этих сил в узлы А и В выполнялся аналогично выражениям (5):

Рв = Р(1); M в = M(Z) + Re х Р(/).

Вектор обобщенных узловых сил КЭВ формировался из проекций векторов сил и моментов, найденных из формулы (7):

Q = (, мА, Pвт, мB)Т.

Для случая единичного обобщенного узлового перемещения с номером г вектор Q представляет собой набор реакций КЭВ на заданное единичное перемещение, т. е. столбец матрицы жесткости. Например, для граничных условий, заданных соотношениями (6), Q дает четвертый столбец матрицы жесткости. Решение каждой линейной краевой задачи позволяет заполнить один столбец матрицы жесткости. Решение двенадцати краевых задач позволило найти все 144 элемента матрицы [К]. При этом контроль осуществляли по выполнению свойства симметрии матрицы жесткости К = Крг.

Построение матрицы масс КЭВ. При решении линейной краевой задачи для системы (1) с единичным вектором Y определялись не только узловые реакции КЭВ, но и распределения векторов перемещений и(5) и поворотов Фг(з) вдоль дуги витка (г — номер единичного перемещения в векторе Y). Объединение линейных перемещений и поворотов {иг (5), Ф г (5)} можно рассматривать как вектор-функцию формы КЭВ с номером г. Подчеркнем, что здесь вектор-функция формы не задавалась, как в традиционном МКЭ, а рассчитывалась интегрированием системы дифференциальных уравнений, поэтому в задачах статики ее можно принять как численно точную.

Элементы матрицы масс [M] определяли как взвешенные скалярные произведения найденных ранее функций формы [10, 19]: I I

тгр = рА|иг • u+ р |фг - JФр (8)

0 0 где г, р — номера функций формы, г = 1, ..., 12, р = 1, ..., 12; р — плотность; I — тензор геометрических моментов инерции сечения.

Как видно из выражения (8), учет инерции вращения осуществляет тензор I. Для круглого сечения тензор инерции вращения может быть представлен в виде, подобном соотношениям (2):

I = ]р +1Х (Е - а).

6x6

¡6x6 1

Рис. 2. Схема сборки матриц [Хг] и ДО^] размером 234x234 из матриц отдельных КЭВ

В отличие от работ [10, 12], в которых рассматривались аналогичные тензоры, здесь элементы тензора J имеют размерность м4 и выражаются через обычные геометрические моменты инерции поперечного сечения.

Сборка матриц жесткости и масс всей пружины и учет закреплений. По найденным матрицам [ЭД и ДО] КЭВ построены результирующие матрицы жесткости и масс ДО£] полной модели цилиндрической пружины. Формирование результирующих матриц пружины из матриц отдельных КЭВ (блоки размерностью 12x12) показано на рис. 2. Как видно из рисунка, диагональные блоки соседних витков размерностью 6x6 суммируются, при этом недиагональные блоки остаются неизменными. Модель пружины содержит 38 витков и 234 степени свободы (по 6 степеней свободы в каждом из 39 узлов).

iНе можете найти то, что вам нужно? Попробуйте сервис подбора литературы.

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

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

веркой. В качестве примера рассмотрена пружина экспериментального стенда [1] со следующими исходными данными:

• модуль упругости материала пружины E = 2 •1011Па;

• коэффициент Пуассона материала пружины ц = 0,3;

• плотность материала р = 8 000 кг/м3;

• радиус цилиндра, описывающего пружину, R = 13,25 •Ю-3 м;

• диаметр проволоки пружины d = 2,6•Ю-3 м;

• угол подъема винтовой линии у = 3,5°;

• количество витков N = 38.

Частоты и формы собственных колебаний пружины находили традиционным для МКЭ способом [20] из однородной системы уравнений

(-p2 [M£] + [K£]) = 0,

где p — искомая круговая частота собственных колебаний; Ys — искомый вектор перемещений (форма колебаний).

Для определения частот и форм к матрице [Ms] [Kе] применялась процедура Eigensystem из программного обеспечения Wolfram Mathematica [17], предназначенная для вычисления собственных значений и векторов матриц.

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

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

Сравнение результатов расчета и эксперимента

Собственная частота Эксперимент Расчет Погрешность, %

f1, Гц 18,3±0,2 18,35/18,34 0,1/0,1

f2, Гц 34,0±0,2 35,82/35,73 5,4/5,1

f3, Гц 39,5±0,2 40,82/40,71 3,3/3,2

f4, Гц 46,6±0,2 46,47/46,42 0,4/0,5

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

0 0,05 0,10 0,15 *ьм

а

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

Сопоставление расчетной и экспериментальной форм колебаний (рис. 3) также показывает их удовлетворительное соответствие.

Выводы

1. Предложена методика расчета частот и форм собственных колебаний «длинных» ци-

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

2. Применение КЭВ позволяет значительно (примерно в 20 раз) сократить размерность конечно-элементной модели пружины и существенно улучшить ее обусловленность.

3. Сравнение результатов расчета, использующих КЭВ, с данными эксперимента показало их хорошее соответствие.

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

Литература

[1] Бадиков Р.Н. Расчетно-экспериментальное исследование напряженно-деформирован-

ного состояния и резонансных режимов вращения винтовых пружин в пружинных механизмах. Дис. ... канд. техн. наук. Москва, 2009. 166 с.

[2] Сорокин Ф.Д., Су Чжоу. Разработка конечного элемента в виде витка цилиндрической

пружины. Нелинейная динамика машин. III Междунар. школа-конференция молодых

ученых (School-NDM 2016), Москва, 12-15 апреля 2016 г., Москва, ИМАШ РАН, 2016, с. 261-265.

[3] Yu A.M., Нао Y. Improved Riccati Transfer Matrix Method for Free Vibration of Non-

Cylindrical Helical Springs Including Warping. Shock and Vibration, 2012, vol. 19, is. 6, pp. 1167-1180.

[4] Yu A.M., Hao Y. Effect of warping on natural frequencies of symmetrical cross-ply laminated

composite non-cylindrical helical springs. International Journal of Mechanical Sciences, 2013, vol. 74, pp. 65-72.

[5] Kacar I., Yildirim V. Natural frequencies of composite cylindrical helical springs under com-

pression. International Conference on Vibration Problems ICOVP, 2011, pp. 119-124.

[6] Yildirim V. On the linearized disturbance dynamic equations for buckling and free vibration

of cylindrical helical coil springs under combined compression and torsion. Meccanica, 2012, vol. 47, is. 4, pp. 1015-1033.

[7] Lee J. Free vibration analysis of cylindrical helical springs by the pseudospectral method.

Journal of Sound and Vibration, 2007, vol. 302, is. 1-2, pp. 185-196.

[8] Becker L.E., Chassie G.G., Cleghorn W.L. On the natural frequencies of helical compression

springs. International Journal of Mechanical Sciences, 2002, vol. 44, is. 4, pp. 825-841.

[9] Hughes T.J.R. The Finite Element Method: Linear statics and dynamic element analysis.

Courier Corporation, 2012. 672 p.

[10] Елисеев В.В. Механика упругих тел. Санкт-Петербург, Изд-во Политехн. ун-та, 2003. 336 с.

[11] Елисеев В.В., Зиновьева Т.В. Механика тонкостенных конструкций. Теория стержней. Санкт-Петербург, Изд-во Политехн. ун-та, 2008. 95 с.

[12] Елисеев В.В., Авксентьев А.И. Модели упругих стержней в динамике гибких роторов. Современное машиностроение. Наука и образование, 2014, № 4, с. 335-343.

[13] Ветюков Ю.М., Елисеев В.В. Моделирование каркасов зданий как пространственных стержневых систем с геометрической физической нелинейностью. Вычислительная механика сплошных сред, 2010, т. 3, № 3, с. 32-45.

[14] Светлицкий В.А. Механика стержней. В 2 т. Т. 1. Статика. Москва, Физматлит, 2009. 408 с.

[15] Феодосьев В.И. Сопротивление материалов. Москва, Изд-во МГТУ им. Н.Э. Баумана, 1999. 592 с.

[16] Валишвили Н.В., Гаврюшин С.С. Сопротивление материалов и конструкций. Москва, Изд-во Юрайт, 2017. 429 с.

[17] Дьяконов В.П. Mathematica 5.1/5.2/6. Программирование и математические вычисления. Москва, ДМК-Пресс, 2008. 574 с.

[18] Попов В.В., Сорокин Ф.Д., Иванников В.В. Разработка конечного элемента гибкого стержня с раздельным хранением накопленных и дополнительных поворотов для моделирования больших перемещений элементов конструкций летательных аппаратов. Труды МАИ, 2017, № 92. URL: http://trudymai.ru/published.php?ID=76832 (дата обращения 14 декабря 2017).

[19] Walter W., Walter D.P. Mechanics of Structures. Variational and Computational. CRC Press, 2002. 912 p.

[20] Ильин М.М., Колесников К.С., Саратов Ю.С. Теория колебаний. Москва, МГТУ им. Н.Э. Баумана, 2002. 272 с.

References

[1] Badikov R.N. Raschetno-eksperimental'noe issledovanie napriazhenno-deformirovannogo

sostoianiia i rezonansnykh rezhimov vrashcheniia vintovykh pruzhin v pruzhinnykh mekhanizmakh. Diss. kand. tekh. nauk [Calculation-experimental study of stress-strain state and the resonant modes of rotation of helical springs in spring mechanisms. Cand. tech. sci. diss.]. Moscow, 2009. 166 p.

[2] Sorokin F.D., Su Chzhou. Razrabotka konechnogo elementa v vide vitka tsilindricheskoi pru-

zhiny [Development of a finite element in the form cylindrical coil springs]. Nelineinaia dinamika mashin. 3 Mezhdunar. shkola-konferentsiia molodykh uchenykh [Nonlinear dy-

namics of machines. 3 International school-conference of young scientists]. 12-15 April

2016, Moscow, IMASh RAN publ., 2016, pp. 261-265.

[3] Yu A.M., Hao Y. Improved Riccati Transfer Matrix Method for Free Vibration of Non-