Математическое моделирование пространственного движения трехзвенного инсектоптера Текст научной статьи по специальности «Механика и машиностроение»

УДК 629.7.014.16

МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ ПРОСТРАНСТВЕННОГО ДВИЖЕНИЯ

ТРЕХЗВЕННОГО ИНСЕКТОПТЕРА

Р.Ю. Поляков, С.А. Бокадаров

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

Ключевые слова: математическое моделирование, пространственное движение, робот.

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

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

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

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

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

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

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

Рассматриваемая система движется в пространстве под действием сил, возникающих при взаимодействии звеньев с воздушной средой. Эти силы приложены в точках А. Кроме того, на трехзвенник действуют силы веса, приложенные в точках С. Для удобства в дальнейшем будем считать, что точки Аi совпадают с С, а крыло представляет собой тонкую однородную недеформируемую пластинку площадью = 2Ьй. Здесь обозначено: 21, й соответственно длина и ширина крыла.

Приближенно сила, возникающая при движении крыла, определяется по формуле

/г= (1/2)С_

где С - безразмерный коэффициент лобового сопротивления, Г - плотность вязкой среды, 5 -эффективная площадь пластинки (крыла), V -скорость движения центра тяжести крыла в воздухе. Анализ этой формулы показывает, что подъемная сила крыла может изменяться в зависимости от скорости движения крыла и его площади. Если управлять этими параметрами по соответствующему

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

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

занимает крайнее положение вверху или внизу. Такой характер движения крыла позволяет изменить площадь крыла и создать асимметрию в законе, определяющем действие подъемной силы и Е3 . Пусть ю=200.. .250 Гц, амплитуда колебаний крыла:

А = 510-3 м.

Приращение площади крыла

8! = 6 10-5 м2.

Плотность воздуха

р = 1,3 кг/ м3

Пусть относительная координата точки приложения подъемной силы крыла движется по закону

= Ас08Ю? .

Тогда скорость движения этой точки определяется по формуле x(t) = - A «sin ct.

Определим силу, возникающую при вертикальном перемещении крыла,

F(t) = а • S(t)• x2 • sign(x), где площадь крыла имеет вид S(t) = S0 + S sin ct.

Средняя за период сила определяется по

формуле

2ж ю

Ff (t d.

2л

Найдем изменение силы для следующих параметров

ю = 200 Гц, 80 = 9 10-5 м2.

На рисунках 1, 2 представлены зависимости

площади 8), скорости х(1) и силы ¥^) для различных значений 8у .

ю = 200 Гц, 80 = 12 10"5 м2, = 0,5

/ 8 '

<4>

res 4 -100

!s^>-0.01 0

<2>

s 2 -2000

/\. S(f)-2( ю/ \ А

А А А

И ° VÍ \ * ■... * 1

\1 \ F-101 V

0.01

<0>

Рис .1. Зависимость площади S(t), скорости x(t) и силы F(t) для Si/ = 0,5

/ S„

ю = 200 Гц, So = 6010-5 м2, Si/ = ОД

0.4

0.2

0.2

0.4

- 3

0

0.015

0.02

res<1>-0.01 о

Sil) ■ 2000

\ h \ /

V V/

0.01

<0>

Рис. 2. Зависимость площади S(t), скорости x(t) и силы F(t) для SiÁ = 0,1

S п

В процессе исследования были получены графические зависимости площади $(/), скорости

х(1:) и силы ¥) для значений равные 0,06;

0

0,1; 0,2; 0,4; 0,5; 0,67; 0,9. Полученные зависимости показывают, что на характер подъемной силы существенно влияет величина . С увеличением

этого параметра до единицы среднеинтегральная сила достигает максимума. А при приближении $ /

к нулю - соответственно стремится к нулю.

На рисунке 3 приведена зависимость максимального усилия от соотношения площадей

V.

0,4 0,6

S1/S0

Рис. 3. Зависимость максимального усилия от соотношения площадей

Рассмотрим силы, действующие на трехзвенник. Пусть на трехзвенник действуют две подъемные силы F1, F3, вызванные взаимодействием крыльев с воздушной средой, приложенные к точкам С1, С3 соответственно, и силы тяжести звеньев т£ , а также силы, возникающие в результате движения воздуха относительно корпуса (возмущающее воздействие).

Силы Fi имеют проекции, параллельные осям Это подъемные силы, удерживающие

объект в воздухе. Векторы этих сил в неподвижной системе координат определим по формулам:

¥(0) = . ¥(1) ¥1 Т12Т20 '

F(0) = TT ■ F(3)

1 3 132120 1 3 '

F <1) = 1 1 =

F <1) = 1 3 =

о о

¥1 ¥3

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

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

ё

ж

Количество движения рассматриваемой механической системы Q определим по формуле

ё=щ >

где

Я = тРс, >

= Z mg+Fi+F3 ■

42 = m2VC

(1)

q 3 = = m3vc3>

vC1 = VC 2 + T20T12 ■ r (D rO1C

VC3 = VC2 + T20T32 ■ r (3) O3C3

+ Т Т Г (!)

+ Т Т г(3)

^ 120132,03С3 . В соответствии с теоремой изменение количества движения равно сумме сил, действующих на 1-ю массу,

dqi dvC =(2)

= m (—С)=t2F + mg ■ dt dt

(2)

Вектор количества движения

рассматриваемой системы с учетом (1), определим по формуле

ё =ЪтРс,. (3)

2

(4> , ЛЛ res ■ 100

2

5x10

0.015

0.02

0

0,2

0,8

2

Теорема об изменении количества движения механической системы в дифференциальной форме имеет вид

dQ

, dv

—— = (m + m + + m1{T10T11rclcl + TwT21rclcl ) +

(4)

dt

dt

+ т (Т Т г <3) + Т Т Г <3)) = V т Р + Т Т ¥(1) + Т Т ¥<3)

-Г ^^*з(Т20Т32'С2С3 — Т 20Т 32' С2С3 ) V '"¡.5 — Т12Т20¥1 — Т 32Т 20¥ 3 ■

Уравнение (4) в проекциях на координатные оси неподвижной системы координат примет вид

(5)

> dv

(m + m2 + m3) + m(T20T12 rC2C1(1) + T20T1 2ГС 2C1(1) ) +

А

+ тз(Т20Т32ГС2С3 — Т20Т32ГС2С3 ) = ^^ тг р — ^12^20^1 — ^32^2(0^3 ■

Если предположить, что масса крыла значительно меньше массы корпуса и привести подъемные силы к корпусу, то в проекциях на оси неподвижной системы координат получим

V = (Бт^т^- соз^соБ^зтб) • V ¥,

Л

Vт С = соБ^т)¥ ,

I т

dV dt

= cospcosd ■ ^F¡ + ^mg

Представленная система дифференциальных уравнений (6) описывает изменение обобщенных координат, определяющих положение центра масс, в пространстве под действием подъемных сил и сил веса.

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

— = У Ме ■ я ^

Так как вектор Ь является подвижным, то производная по времени определится следующим образом:

= ^ + ) = (7)

Кинетический момент элементов

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

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

скорости (2) крыльев (1 = 1, 3) представим в виде векторной суммы

П(2)г = Ю (2) + ю<2) + ю(2),

где

а,

(2)

ä

ä

ä

вектор угловой скорости

корпуса, определенный через проекции на оси системы координат C2X2Y2Z2,

=

ä.

0

- вектор угловой скорости

крыльев 1 = 1, 3, определенный через проекции на оси системы координат О1Х1У121,

о . = а - проекции векторов угловых скоростей крыльев на оси 02Х2, О4Х4. 0

С) =

0

ä

- вектор угловой скорости

крыла в системе координат OjXjYjZj.

Определим векторы угловых скоростей во второй системе координат.

— (2) ^ —(i) „ С. = 1 а ■ с. - вектор угловой скорости

крыльев в проекциях на оси системы координат C2X2Y2Z2.

Здесь обозначено:

'Т"' _ ГТ7 ф 2

1 i2 = 1i 2 1 i 2 '

где

cos a 0 sin a

0 1 0

- sin a 0 cos a

от системы координат OjXjYjZj к системе координат CX2Y2Z2 при повороте на угол а.

1 0 0 1 = cos^ sm^2 0 - матрица перехода

- sin a cos фп 0

от системы координат OjXjYjZj к системе координат CX2Y2Z2 при повороте на угол фа . Отсюда

гра _

2 =

- матрица перехода

cosa 0 sina 0 ä sina

ä(2) = 0 1 0 0 = 0

- sina 0 cosa ä cos a

wY X 2 cozl sina а X 2 + äzi sina

aY 2 + 0 = aY

äZ 2 ®Z2 cosa äZ 2 + czi cosa

Тогда в частном случае при а = 0 ю^ = 0

абсолютная угловая скорость вращения крыльев имеет вид, что соответствует формуле

П(2) =

Теперь найдем момент количества движения рассматриваемой механической системы

Ь = Ьс 2 Ь ,

где ЬС2 = Iс2юС2 - кинетический момент

корпуса относительно центра масс инсектоптера;

Д = /. ^ - кинетический момент г-го

крыла относительно центра масс инсектоптера. Тензор инерции корпуса 1С2 равен

0

JX 2 J С 2 0 0

2 = 0 ,/r 2 J С 2 0

0 0 ^2 J С 2

где jX 2 JC2 , ./7 2 J С 2 tZ 2 , JC2

- осевые моменты

инерции корпуса.

Если упрощенно предположить, что корпус является параллелепипедом с размерами l1, l2, l3 по соответствующим осям, то осевые моменты инерции можно определить по формулам:

JXc2 = 1 + £) > JYcl = 1 m(& + /32) >

J,

J_ 12

С2 = —m(l + l2).

Для вычисления тензора инерции крыльев для системы координат С2Х2У222 построим вспомогательные системы координат 01Х'1У'12'1 параллельные системе координат С2Х2У222.

Моменты инерции крыльев в системе координат 01Х1У121 определим по формулам:

^ = +12),

3

JO = m(li +132) .

JÖt =- m(l2 +,

Направление осей OjXjYjZj (i = 1, 3) в системе координат C2X2Y2Z2 будет определяться таблицей направляющих косинусов, образующих матрицу:

аи а12 а1Ъ

= Му (а) •

Вид Му (а) определяется типом вращения

одной системы координат относительно другой. Поскольку в исследуемой механической системе оси вращения крыльев OZi могут осуществлять поворот на углы а только относительно OYi, получим искомую матрицу в виде

(cosa 0 sina^

'aii ai2 ai3N

a2i a22 a23

va3i a32 a33 ,

где

М у (р) =

0 1 0 - sina 0 cosa

- матрица

поворота осей крыла относительно оси СУ1.

Определим теперь осевые и центробежные моменты инерции несущих винтов относительно системы координат СХ^^:

ГА1 = т(у2а + 22а) + + с^ + аЩ,

I* = т^2а + хС,)+ а^Ц + с^ + а^А,, !а, = т{х2с, + у2С1)+ + ^ + аЦА

= тХаУа + ai ia2 1 (i

Iy = тУсггсг +a2 iasi (i

= юга*« +a ia3i(l

TXl

- iL

Xl

- IAl

+ ai 2a22(IA4

(I (I

)+ai 3ai 2(i

+ a22a32(iA

Iyl) Iyl)

IAi)'

Iyl)

I Ai),

где

IAi = 0, IAyi = ^ m(3r2 + h2)

1

/Zl 2

лг = mr - моменты инерции несущих винтов

относительно осей системы координат OXYZ m - масса i-го крыла,

Аа = 0, yCi = l, zc¿ = 0 - координаты i-го крыла в системе координат CXYZ,

an = cos Р, a12 = 0, a13 = sin Р, a21 = 0, «22 = 1, a23 = 0, ««i=- sin Р, а32 = 0, a33 = cos Р.

Отсюда получим: ГА1 = ml2 + sin2 a • IZI,

Iy = Iy

1 Al 1 Al ,

Гм = ml2 + cos2 a • I^.,

IAy = 0

1 Ai' 0 ,

Iyz = 0

I Al 0 ,

^ = cosa(" sin a)I% . Таким образом, тензор инерции i-го крыла Ii примет вид

0 - cos р(- sin р)1 mr2

ml2 + sin2 P- — mr2 2

0 — m(3r2 + h2) 0

12 V У

- cos р(- sin р)"1 mr2 0 ml2 + cos2 р-mr2

Кинетический момент корпуса относительно центра масс

JX 2 J С 2 0 0 ^X2 JC2

Lc = Icac = 0 JY 2 2 0 = rY 2 •с 2-щ7 2

0 0 •Z 2 J С 2 щ Z 2 rZ 2 •с2 -WZ 2

Кинетический момент i-го крыла относительно центра масс

ml2 + sin2 P- - mr2 2

- cos р(- sin р) i mr2

—m(3r2 + h2) i2 v '

ml2 + cos2 P -1 mr2 2

yi

щ cos a

(ml2 + sin2 P-imr2)-(щ + щ sinP)-imr2 cos2 р(-sinр)-щ

— m(3r2 + h2 )-i2

y

-1 mr2 cosр(- sinР)(^ + ^ sinР) + (ml2 + cos2 Р-1 mr2)cos Р-.

Теперь определим кинетический момент трехзвенника в целом

(ml + sin2 a •1 mr2) • (шх + шг; sina) -1 mr2 cos2 a(- sina)-ü

jX 2 JC 2

J'd +S

JZ 2

2 x « 2 —m(3r2 + h2 )-<

n v '

imr2 cosa(- sina)(^ + ^ sina) + (ml2 + cos2 a -1mr2)cosa - шг

Jc&Xl + (ml2 + sin2 a-1 mr2) - (2^ + (^ + ®г4 )sina) -1 mr2 cos2 a(- sina)(^ + ^) J£+ ^ m(3r2 + h2)- (®y2 +®y4)

Jr, — mr cosa(-sina)(2^

^y2 + ^y4)

(-sina)(2^ + (^ + ®z4)sinP) + (ml2 + cos2 a-imr2)cosa- (^ + )

cos р(- sin р)—mr2

0

(°x¡ +Щ sina

Li = I"i

0

0

0

L

L

Теорема об изменении кинетического момента механической системы при движении относительно подвижной точки С2

— = 7 Щ. <и 7 с

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

dL dt

dL dt

dL dt

Jx<

Jy<y +

J z< z

JX<x

Jy<y +

JZ<z

JX<x

Jy<y +

JZ<z

i <x

JXo„

J

<y

k

<z

Jz<

J2aya2 - Jy<y<2 J<<z - Jz<x<z

Jy« - Jx

<x<z

<x<y

Jz - J

(Jx - J Jy - J

M

y 2

(8)

(9)

M

z

x

Подставляя (9) в (7) найдем об изменении кинетического момента в проекциях на

соответствующие уравнения, отражающие теорему связанную систему координат:

' Гюх +ЮуЮг (Г - ^ )=

<Гюу +юхюг (Г - Jz )=Му2, (10)

Jzюz -ЮхЮу(/у - Jx)=М22.

Формулы (6) и (10) позволяют построить систему дифференциальных уравнений,

описывающих движение инсектоптера.

Разработаная математическая модель позволяет описывать пространственное движение

Библиография

1. Красковский Е.Я. Расчет и конструирование механизмов приборов и вычислительных систем: учеб. пособие для приборостроит. спец. вузов / под ред. Ю.А. Дружинина. - 2-е изд., перераб. и доп. - М.: Высш. шк., 1991.

2. Поляков Р.Ю. Изучение движения взлета летающего робота с машущим крылом / Р.Ю. Поляков, С.В. Ефимов, Р.И. Праслов // Вестник Воронежского института МВД России. - 2014. - № 3 - С.90-97.

3. Поляков Р.Ю. Исследование управляемого синхронного движения летающего многозвенного робота /Р.Ю. Поляков, С.В. Ефимов, Н.В. Мозговой

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

References

1. Kraskovskij E.YA. Raschet i konstruirovanie mekhanizmov priborov i vychislitel'nyh sistem: ucheb. posobie dlya priborostroit. spec. vuzov / pod red. YU.A. Druzhinina. - 2-e izd., pererab. i dop. - M.: Vyssh. shk., 1991.

2. Polyakov R.YU. Izuchenie dvizheniya vzleta letayushchego robota s mashushchim krylom /R.YU. Polyakov, S.V. Efimov, R.I. Praslov // Vestnik Voronezhskogo instituta MVD Rossii. - 2014. - № 3 -S.90-97.

3. Polyakov R.YU. Issledovanie upravlyaemogo sinhronnogo dvizheniya letayushchego mnogozvennogo robota / R. YU. Polyakov, S. V. Efimov, N. V. Mozgovoj // EHlektrotekhnicheskie

// Электротехнические комплексы и системы управления, Воронежский инновационно-технологический центр. Воронеж: Кварта, 2014. -№ 3 - С 28-33.

4. Яцун С.Ф. Исследование движения многозвенного робота, перемещающегося прыжками и планированием / С.Ф. Яцун, Л.Ю. Волкова, А.В. Ворочаев // Справочник. Инженерный журнал с приложением. - 2014. - № 4. С. 12-17.

kompleksy i sistemy upravleniya, Voronezhskij innovacionno-tekhnologicheskij centr. Voronezh: Kvarta, 2014. - № 3 - S 28-33. 4. YAcun S.F. Issledovanie dvizheniya mnogozvennogo robota, peremeshchayushchegosya pryzhkami i planirovaniem / S.F. YAcun, L.YU. Volkova, A. V. Vorochaev // Spravochnik. Inzhenernyj zhurnal s prilozheniem. - 2014. - № 4. S. 12-17.

MATHEMATICAL MODELING OF SPATIAL MOVEMENT OF THREE-TIER

INSERTAFTER

This article discusses the mathematical model of the movement of insectopia based on the three-tier model of Electromechanical system with two oscillatory motion of the outer links, leading to the formation of vibrational effects, as by inducing a lifting force and pulling force realized by the use the effect of "asymmetry" wing shape and speed. Particular attention is paid to the mathematical description of the properties of electric drives, the kinematics of the rotation of the body, algorithms for the formation of control voltages, gyroscopic moments of rotating rotors of electric motors and external periodic perturbations.

Key words: mathematical modeling, spatial motion, robot.

Поляков Роман Юрьевич,

преподаватель кафедры технологических процессов в машиностроении и агроинженерии,

ФГБОУ ВО «Елецкий государственный университет им. И.А. Бунина», e-mail: polyakov_gps@mail.ru, Россия, г. Елец, Polyakov R. Y.,

lecturer of the Department of technological processes in engineering and Agroengineering, Bunin Yelets State University, Russia, Yelets.

Бокадаров Станислав Александрович,

кандидат технических наук,

научный сотрудник организационно-научного и редакционного отдела,

ФКОУ ВО Воронежский институт ФСИНРоссии,

e-mail: bokadarov.stas@inbox.ru,

Россия, г. Воронеж,

Bokadarov S.A.,

candidate of engineering science,

researcher of organizational, scientific and editorial Department, Voronezh institute of the Russian Federal penitentiary service, Russia, Voronezh.

© Поляков Р.Ю., Бокадаров С.А., 2018