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

Том XXXIII

УЧЕНЫЕ ЗАПИСКИ ЦАГИ 20 0 2

№ 1—2

УДК 629.7.015.017.2

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

А. О. БЕЛЯКОВ

Представлено компьютерное исследование нового метода точного (+0,1%) определения моментов инерции крупногабаритных тел, разработанного в ЦАГИ. Метод предполагает свободные колебания твердого тела по трем степеням свободы. Предложен способ моделирования процесса колебаний твердого тела с несколькими степенями свободы. Определены условия применимости метода Прони нахождения собственных частот сигнала, полученного в результате моделирования. Приведены результаты моделирования колебаний конкретной системы при различных положениях центра масс.

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

ср + -|ф=0,

Р

где г — расстояние от оси закрепления до центра масс; g — величина ускорения свободного падения; р — радиус инерции, который связан с моментом инерции I соотношением I = Мр , где М — масса тела.

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

2 2 2 Рц .м =Р - Г .

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

В. В. Богдановым был разработан новый метод измерения моментов инерции по свободным колебаниям тела по трем степеням свободы. На основе этого метода в НИО-7 ЦАГИ был

Рис. 1. Схема измерительного стенда

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

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

1. Численное моделирование процесса колебаний. Для исследования процесса колебаний написана программа, моделирующая движение системы путем численного решения уравнений динамики:

= 1 *,

Л

2

(1.1)

(1.2)

т ъ» ^2ф

где I — тензор инерции, ф — вектор, компонентами которого являются углы Эйлера, —— —

вектор углового ускорения, у — радиус-вектор центра масс,

&2 у & 2

Л2

вектор ускорения центра

масс, г — радиус-векторы точек приложения сил ^. Так как силы ^, за исключением силы тяжести, являются силами упругости, то движение будет колебательным. Система, описываемая уравнениями (1.1) и (1.2), имеет, вообще говоря, 6 степеней свободы, но, приложив силу

параллельно оси У, мы можем возбудить колебания только по трем степеням (см. рис. 1). Кроме того,

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

Опишем процесс моделирования колебаний. Он представляет собой численное решение дифференциальных уравнений второго порядка (1.1), (1.2) при помощи модифицированного метода Эйлера [3], [4] с однократным пересчетом. Ниже приведена последовательность операций на одном шаге численного решения длительностью А?. Рассматриваются две системы отсчета: система, связанная с телом с началом координат в центре масс тела и осями вдоль главных осей инерции, и неподвижная система отсчета, показанная на рис. 1. В начале процесса пружины недеформированы и соответствующие оси обеих систем сонаправлены, а векторы ускорения центра масс тела а, углового ускорения б, поступательной V и угловой ю скоростей нулевые.

Последовательность операций на «-ом шаге моделирования:

1) производится расчет положения подвижной системы отсчета по прошествии вре-

А/

мени —:

2

Аt А/12

Уп+1/2 = Уп + + ап~,

А / М2

ф п+1/2 =ф п +ю пу+б п—;

2) находится матрица поворота Б на угол фп+1/2 -фп системы отсчета, связанной с телом, относительно неподвижной системы;

3) определяются новые радиус-векторы от центра масс к точкам крепления пружин в неподвижной системе отсчета:

Г1, п+1/2 = ^ ' Г 1,п + у и+1/2’ , ^ = 1, 2, 3, 4;

4) из матрицы Б и геометрических соотношений находятся векторы сил ¥ г-п+1/2,

I = 1, 2, 3, 4;

5) определяются векторы ускорения центра масс в неподвижной системе отсчета:

_ 1

а п+1/2 ¥ 1, п+1/2 + 3;

1 =1

где g — вектор ускорения свободного падения;

6) находится тензор инерции в неподвижной системе отсчета:

+1/2 = $Т !п $;

7) определяется вектор углового ускорения £ п + 12 с помощью решения системы линейных уравнений:

4

!п +1/2£п +12= 2¥ 1+1/2 х г 1 +1/2;

1=1

8) находятся положение и скорости подвижной системы отсчета по прошествии всего интервала Аt, используя ускорения, вычисленные на его середине:

Уп + 1 = Уп + VnAt +

an + 1/2 At 2

v n+1 = vn + a n + 1/2^

9) вычисляется матрица поворота Б на угол ф п +1 - ф п +1/ 2;

10) затем производятся операции пп. 3—7) для определения векторов ап +1 и £ п +1 и т. д.

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

v0 . , л . ГС

аналитическим решением y(t) =—sin(®t) + y0 cos ), где ю = . —, с — жесткость пружины,

ш \М

можно заключить, что точность данного численного метода лежит между вторым и третьим

порядком мелкости разбиения.

Благодаря относительной простоте метод Эйлера с пересчетом показывает, что при

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

характере движения тела. Для аппроксимации движения на каждом шаге алгоритма

используются только уравнения равноускоренного движения и второй закон Ньютона, что

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

натурному эксперименту на ранних этапах конструирования подобных измерительных стендов.

2. Определение частот, декрементов затухания и амплитуд сигнала методом Прони.

Если число мод колебаний в сигнале известно заранее, то их частоты и декременты затухания

можно определить с большой точностью, используя метод Прони [2]. Вначале коротко опишем

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

p a

x [n] = Ao +^ -k [exp (rnknT ~bknT + % ) + exp (-mknT-bknT - % )] , (2.1)

k=1 2

где P — число мод колебаний в сигнале, —к — амплитуда к-й моды, 5k — декремент затухания к-й моды, Шк — частота к-й моды, фк — фаза к-й моды, T — период дискретизации сигнала.

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

2 P+1

^ акХ[n - к] = 0, (2.2)

к=0

где 2P +1 < n < N — 1, N — число элементов в массиве х. Чтобы найти коэффициенты ак, нужно решить систему из 2P +1 линейных действительных уравнений, положив ао = 1. Эта система получается из условия:

2

min

ak (k=0,2 P+1)

N-1 f 2 P+1 V

Z Z akx[n - k]

IV k=0

n=2P+1

Далее находятся корни следующего полинома степени 2P +1:

2 P+1

F(z)=Z akz2P+1-k. (2.3)

k=0

Один из его корней является вещественным, а остальные составляют комплексно-сопряженные пары:

zk = exp (mknT -5knT), z-k = exp (-i&knT-bknT), к = 1, ..., P. (2.4)

Таким образом, из корней (2.4) полинома (2.3) мы находим частоты и декременты затухания колебаний.

Для определения значений hk = A- exp (iфk) нужно решить систему из 2P + 1 линейных

уравнений с комплексными коэффициентами.

При применении метода обнаружено, что он дает очень точные (~ 0,001%, N = 100) значения частот и декрементов затухания слабо зашумленного сигнала (белый шум 1%), при

T|max е [2, «10), (2.5)

где ®max — максимальная частота составляющих сигнала, причем точность повышается, при

приближении значения T ш max/(2л) к левому краю интервала в соотношении (2.5). Из теоремы

Котельникова [5] также следует, что левая граница интервала не может быть меньше 2. В то же время для определения амплитуд мод колебаний частоту разбиения сигнала по времени лучше выбирать так, чтобы на период моды 2л/шх, амплитуду которой мы хотим определить, приходилось бы примерно столько же точек, сколько периодов этой моды укладывается в суммарное время измерения N • T, т. е.

T . <2-6>

Отсюда следует, что при больших N соотношение (2.5) может не выполняться. В таком случае для определения частот нужно рассматривать не все точки сигнала, а каждую m-ю

m , (2.7)

1 ^nax

где T взято из выражения (2.6). Тогда соотношение (2.5) будет выполнено, так как вместо T в (2.5) будет Tm.

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

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

Вычисления проводились для тела диаметром 1,74 м, подвешенного на пружинах длиной 0,5 м (см. рис. 1). Расстояние между пружинами по оси X составляло 4 м, жесткость пружин равнялась 740 250 Н/м, а масса тела — 15 190 кг. Главные моменты инерции равнялись 3703,

43 920, 43 920 кг • м2. Колебания возбуждались в точке крепления одной из пружин приложением силы

0,4

Рис. 2. Распределение энергии по модам колебаний (Дж)

Рис. 3. Амплитуды колебаний по модам и номерам креплений пружин

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

Г 4

| Е Л?,.

2 г=1

где С — жесткость пружины, Лк і — амплитуда к-й моды колебаний на 7-й пружине, к = 1, 2, 3, 4.

При различных положениях центра масс тела в плоскости ОХХ было посчитано 999 точек (рис. 4), что заняло « 20 часов времени работы компьютера (Pentium 1, 160 МГц). Из полученной поверхности минимумов энергий в модах колебаний (рис. 4) видно, что при возбуждении колебаний силовым воздействием в точке крепления одной из пружин три моды колебаний хорошо определяются, если смещения центра масс относительно центра жесткости не превышают ~ 20% от расстояния между пружинами по соответствующей оси. На рис. 5

изображена зависимость двух частот системы от смещения центра масс относительно центра жесткости по оси X. Так как смещение по оси Z равно нулю, то частота вращательных колебаний вокруг оси X остается постоянной и равна 24,59 рад/с. На рис. 5 видно, что с увеличением смещения частота ю колебаний системы вдоль вертикальной оси монотонно уменьшается, а

частота вращательных колебаний вокруг оси Z монотонно увеличивается; это можно объяснить физически. Полученные

из этих частот моменты инерции совпадают с изначально заданными моментами инерции тела с точностью ~ 0,1%.

Рис. 4. Минимальные энергии мод колебаний (Дж)

Смещение центра масс относительно центра жесткости по оси X(м)

Рис. 5. Частоты по модам колебаний (рад/с)

Результаты численных опытов показывают, что измерительный стенд представленной конструкции имеет достаточно широкие пределы по смещению центра масс тела относительно

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

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

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

ЛИТЕРАТУРА

1. Маркеев А. П. Теоретическая механика.— М.: Наука.— 1990.

2. Марпл - мл. С. Л. Цифровой спектральный анализ и его приложения. — М.: Мир,

1990.

3. Бахвалов Н. С., Жидков Н. П., Кобельков Г. М. Численные методы.—

М.: Наука.— 1978.

4. Годунов С. К., Рябенький В. С. Разностные схемы. Введение в теорию.—

М.: Наука.— 1977.

5. Темников Ф. Е. Афонин В. А. Дмитриев В. И. Теоретические основы информационной техники.— М.: Энергия.— 1971.

Рукопись поступила 16/Ш 2001 г.