Имитационное моделирование спуска осесимметричного авторотирующего тела в квазистатической среде Текст научной статьи по специальности «Физика»

Беляков Д.В.

МАТИ-РГТУ им. К.Э. Циолковского, Москва, Оршанская ул., д.3, к.т.н., доцент кафедры Прикладная математика, информационные технологии и электротехника (495)-141-95-

57, pm@mati.ru

Имитационное моделирование спуска осесимметричного авторотирующего тела в квазистатической среде

КЛЮЧЕВЫЕ СЛОВА

Тело сложной конфигурации, стационарный режим, осреднение, имитационное моделирование.

АННОТАЦИЯ

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

Введение

Быстро вращающиеся тела обладают целым рядом важных свойств, которые использовались человеком в очень давние времена. Эти свойства использовались при охоте с помощью бумеранга, при ведении боевых действий, в детских игрушках-волчках. Классические исследования по динамике абсолютно твердого тела Эйлера и Лагранжа дали обоснование некоторых свойств вращающегося твердого тела. Всем известно, что быстро вращающийся волчок обладает свойством сохранения направления своей оси, также как и ротор гироскопа. Такие же свойства имеет авторотация, т.е. свободное вращение воздушного винта под действием набегающего потока воздуха. Полет вертолета, различные ветродвигатели и многие другие вещи связаны с авторотацией. В 1920 году Дарье предложил идею ветродвигателя с вертикальной осью, который вращается в горизонтальной плоскости. Преимущество такой конструкции перед другими турбинами состоит в том, что она не требует вмешательства в свою работу при смене направления ветра. В то же время, у ветротурбины Дарье момент вращения намного больше чем у обычной крыльчатой ветротурбины (винта). Но если обычный несущий винт известен как простейшая система спуска в режиме авторотации, то неизвестно как конструкцию ветряка Дарье использовать в качестве «парашюта».

Постановка задачи

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

стержню (Рис. 1). Будем считать, что движение происходит в вертикальной плоскости. Выбор в качестве предмета исследования достаточно простого механического объекта-элемента ветротурбины Дарье с горизонтальной осью оправдан с точки зрения интереса к его использованию в качестве системы спуска.

При создании модели воздействия среды на тело используется гипотеза о квазистационарном обтекании пластинок средой [10], [11], [12], [13]. Согласно этой гипотезе сила воздействия среды на каждую пластинку характеризуется скоростью некоторой ее точки, которая называется центром давления. Необходимо отметить, что подобный способ описания аэродинамических сил, действующих на тело, не является единственным. Помимо упомянутой используются вихревая [14] и многотрубчатая [16] модели взаимодействия пластинок со средой. Активно исследуется влияние неоднородности потока на аэродинамические нагрузки [15] и эффект срыва потока. Выбранная квазистатическая модель обладает рядом преимуществ, таких, как простота, наглядность, хорошая согласованность с экспериментом. В рассматриваемой модели предполагается, что поперечные размеры пластинок намного меньше длины стержня. При выполнении этого условия, центры давления пластинок точки А и В можно

I

Рис. 1

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

Составим, уравнения движения рассматриваемого тела в системе координат, связанной с осями Кенига. При этом суммарные составляющие аэродинамических сил, действующих на каждую пластинку, удобнее представить в виде суммы касательных сил ТА, Тв, направленных против пластинок и нормальных сил ЫА, Ыв, ортогональных им (Рис. 1).

Величины касательной и нормальной сил равны:

| TA |= r(a + ó)VA2 = 0.5pscT(a + ó)VA2,

A ^ i ^ /'A

т2 п г___ г п . j\тл2

Na |= n(a+8)V¡ = 0.5pacn(a + S)V]

A у A

22

| Тв |= т(Р + 5УВ = 0.5рает(Р + 5УВ, | Ыв |= п(Ь + 5)УВ = 0.5роеп(Ь + 5)УВ

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

т(а+5) = s(а + 5)cos(а + 5) - р(а + 5) sin(а + 5), п(а + 5) = s(а + 5^т(а +5) + р(а + 5)cos(а +5)

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

Тогда теорема о движении центра масс в проекциях на оси Кенига и теорема об изменении кинетического момента будут иметь вид:

mVx = (t(a + 5)V¡ -t(fi + 5)VB2)sin(y + 5) + (n(fi + 5)V¡ -n(a + 5)VA2)cos(y + 5) + mg mVy = (t(b + 5)VB2 -t(a + 5)VA2)cos(y + 5) + (n(b + 5)V¡ -n(a + 5)VA2)sin(y + 5)

Jó) =-rcos5(t(a + 5)Va2 + t(b + 5)VB2) + rsin5(n(fi + 5)V¡ -n(a+5)V¡) y) = ó

(1)

(2)

Кинематические соотношения, связывающие VA,VB,a, b с Vx,Vy,¥,y) , имеют вид:

VA sin a = (Vy sin y + Vx cos y)

VB sin b = -(Vy sin y + Vx cos y),

VA cos a = ry - (Vx sin y- Vy cos y)

VB cos b = ry) + (Vx sin y - Vy cos y)

Аэродинамические функции T(ff), J<-ff) связаны

соотношениями:

t(a + 5) = s(a + 5) cos(a +5) - p(a + 5) sin(a + 5) n(a +5) = s(a + 5)sin(a +5) + p(a + 5)cos(a +5) После того, как мы проинтегрируем систему уравнений (1) - (2),

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

x = К,У = V. (3)

Таким образом, построена математическая модель движения тела, представляющая замкнутую систему уравнений (1)-(3).

Простейшие установившиеся движения

Будем искать простейшие установившиеся режимы движения, при которых у тела отсутствует вращение. При этом тело будет совершать поступательное движение, т.е. w □ 0. Правая часть уравнения вращения системы (1) обращается в нуль тождественно, т.к. Jw = 0.

Действительно, аэродинамические функции s(a) p(a) имеют период л, а углы атаки при поступательном движении связаны соотношением: ß=л+a.

Проводя простейшие преобразования, получим:

Jw = -V2 r cos S (t(a + S) + t(ß+S)) + V 2r sin 5(n(a + S) + n(ß + S)) = = -V2r cos S(t(a + S) - t(a + S)) + V2r sin S(n(a + S) - n(a + S)) □ 0

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

Таким образом, любое стационарное решение решение вида y(t) □ const = y* является решением уравнения вращения и является установившимся движением.

Введем для удобства угол отклонения вектора скорости от вертикали и обозначим его У. Назовем этот угол углом планирования.

Из геометрических соображений:

у = у+а-Л (4)

Для определения стационарного значения | V | запишем уравнения движения центра масс:

2t (a +S)V2 sin(y + S) - 2n(a+ S)V 2cos(y + S) + mg = 0

2t (a + S)V 2cos(y +S) - 2n(a + S)V2 sin(y + S) = 0 (5)

После несложных преобразований перепишем второе уравнение системы (4) в виде уравнения:

ctg (y+a) = k (a+S) (6)

k ( ) n(a)

где k (a) = -—.

t(a)

Теперь мы можем задать произвольное стационарное решение уравнения вращения y(t) □y = const и при помощи решения уравнения (6) и использования соотношения (4) определить a . Таким образом мы полностью определили направление скорости V.

Рис. 2

Из первого уравнения (5) определим величину скорости:

=J 8 (7)

В работе [3] дана графическая интерпретация решений уравнения (6) и приведены возможные стационарные режимы. Зависимость g = g(w) угла планирования g от угла установки тела y имеет вид (рис.2.).

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

Определим режим движения при высокой угловой скорости. Будем считать, что ™ » v и тело авторотирует. При а -» О и Р -» 0 линеаризуем кинематические соотношения:

sillff'^ff sin , COSff^l И СО50™1 VAa = (l^cosijj + Vy sin = —(I^coslJj + Vy sin ф)

VA = r(á — (l^ sin ф — Vyc<os.ty}VB = t(ií + sin iJj — Vy созф) Аэродинамические функции

и «(«) линеаризуем при а -»О и

>6 —> 0 :

s(a) = ^p<J(cte +с[5а),р{а) = ^poic^

В уравнении движения и уравнениях движения центра масс сделаем

переход от переменных Ул, ^в, а и Р к переменным Ух, ^у, ф и после чего проведем их осреднение по быстрой переменной у. Полученная система будет иметь следующее стационарное решение:

_ -¡К ¡Ср-Схе „ __УЩ_ „ т%(с[6-Ъсуе)8

^ ^ -рогщ^ + Зс,,) -

-рогщ^ + Зс^) -рогщ^+Зс^У

(8)

Это установившееся движение носит название режима авторотации.

Скорость вертикального снижения в таком режиме существенно меньше угловой скорости.

Сравнение установившейся скорости режима авторотации со скоростями других режимов

Проведем сравнение стационарного значения вертикальной составляющей скорости на режиме авторотации со стационарными значениями вертикальных составляющих скоростей на других установившихся режимах, при которых тело движется поступательно. Установившийся угол планирования в режиме авторотации у0 пропорционален величине угла перекоса пластинок 5 . Если он достаточно мал, то вертикальная составляющая скорости в первом приближении будет иметь такой же вид, как в формулах (3.1.5), определяющих параметры режима авторотации.

2 I 2с.

гО

I Сх0

2 1

С другой стороны, для минимального значения вертикальной составляющей скорости при поступательном движении в § 2 Главы 1 при помощи первого приближения была получена оценка (2.1.14). Таким образом, достаточно сравнить:

1 2СГ0 „ |

И

Для тела, имеющего прямоугольные пластинки с удлинением (1 = 8), имеем:

^ = 0.01,с"0 Я 2.62423,сх2 = = 1.312115,^ = 4.58.

1 = 0.0379 и |-^-^ = 0.2498

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

Имитационное моделирование множества стационарных

режимов.

В работе [3] было проведено аналитическое исследование поступательного движения тела. В результате поиска графических решений уравнения (6) при различных значениях угла ориентации тела y качественным образом была изображена зависимость g = g(y) (рис. 2). Но графический метод исследования позволяет найти лишь некоторые свойства решения и лишь приблизительно оценить их. Для более точного поиска корней уравнения (6) написана программа, с помощью которой ищется численное решение этого уравнения при различных значениях y и 8 = 0.

Задача определения точек пересечения графиков y = ctg(y + a) и

y = k (a) сводится к поиску численных решений системы уравнений:

jctg (y+a) - y = 0

(a) - y = 0

При решении этой системы необходимо учесть, что функция

k t ч p(a)

качества y =k (a) = Sa) для плоских прямоугольных пластинок, как

правило, представляет собой экспериментальные данные. Например, для

прямоугольных пластинок с удлинением 1 = 8, такие данные на интервале

pp

[- "2>"2"J, имеются в [13].

Для решения систем нелинейных уравнений в среде программирования MATLAB используется процедура fsolve, входящая в пакет расширения Optimization Toolbox. Эта процедура позволяет решить системы нелинейных уравнений уравнений вида: F (X) = 0 методом наименьших квадратов. Необходимо подготовить файл-функцию для этой системы уравнений. Преобразованная система уравнений для составления файл-функции будет иметь вид:

F (х) = JfXl'Х2) = ctg (х(!) + У) - х(2)

\f2(xi, x2) = spline(alfa,k,x(1)) - x(2) Значение функции качества в каждой точке итерационного процесса при поиске численных решений будем приближать кубическим сплайном.

После этого составим в редакторе MATLAB файл-функцию следующего вида:

function f=myfyn3(x) f(1)=cot(x( 1)+ksi)-x(2) f(2)=spline(alfa,k,x(1))-x(2) Далее вводится начальная точка итерационного процесса и происходит вызов процедуры fsolve: a=spline(alfa,k, 0.01 )

x=fsolve(@myfyn3,[ 0.01 a ],optimset('Display','off)) После завершения итерационного процесса найденное численное

решение выводится с точностью, близкой к машинной. В качестве результата моделирования выводятся графики кривых 7 = у(у) и V = Ух(у)

Рис. 3

Таким образом, задав положение тела при помощи угла у и установочный угол пластинок при помощи угла 5, мы можем численно найти стационарные значения V, у0 при поступательном движении тел

Заключение

Таким образом, в работе:

1. Создана математическая модель движения тела.

2. Получены следующие результаты:

- найдено множество неизолированных установившихся режимов, при которых у тела отсутствует вращение;

- найден стационарный режим авторотации. Показана зависимость режима авторотации от аэродинамических характеристик;

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

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

Литература

1. Локшин Б.Я., Привалов В.А., Самсонов В.А. «Введение в задачу о движении точки и тела в сопротивляющейся среде». Издательство Московского университета. 1992.

2. Беляков Д.В., Самсонов В.А., Филиппов В.В. «Исследование движения несимметричного тела в сопротивляющейся среде» // Издательство «МЭИ», журнал «Вестник МЭИ», выпуск № 4 2006 г., стр. 5-10.

3. Беляков Д.В. «Исследование и особенности математической модели движения несимметричного авторотирующего тела в квазистатической среде» // Издательство «Новые технологии», журнал "Мехатроника, Автоматизация, Управление". Выпуск № 11. 2007 г., стр. 20-24.

4. Самсонов В.А., Беляков Д.В., Чебурахин И.Ф. «Вертикальное снижение тяжелого симметричного авторотирующего тела» в сопротивляющейся среде» // Издательство «МАТИ»-РГТУ, сборник «Научные Труды МАТИ», выпуск 9 (81). Москва. 2005. г., стр. 145150.

5. Самсонов В.А., Беляков Д.В. «Математическое моделирование движения симметричного авторотирующего тела, раскрученного до высокой угловой скорости, в воздушной среде». Издательство «МАТИ» -РГТУ сборник «Научные Труды МАТИ» выпуск 10 (82). Москва. Изд-во МАТИ-РГТУ 2006 г., стр. 196-200.

6. Табачников В.Г. «Стационарные характеристики крыльев на малых скоростях во всем диапазоне углов атаки». Труды ЦАГИ 1974 г. выпуск 1621.

7. Журавлев В.Ф., Климов Д.М. «Прикладные методы в теории колебаний». Издательство «Наука» 1988 г.

8. Волосов В.М., Моргунов Б.М. «Метод осреднения в теории нелинейных колебательных систем». Издательство Московского университета. 1971.

9. Боголюбов Н.Н., Митропольский Ю.А. «Асимптотические методы в теории нелинейных колебаний». Издательство «Наука». 1974.

10. Малкин И.Г. «Теория устойчивости движения». Издательство «Наука» 1966 г.

11. Тихонов А.Н., Васильева А.Б., Свешников А.Г. «Дифференциальные уравнения». Издательство «Наука» 1985 г.

12. Бахвалов Н.С., Жидков Н.П. Кобельков Г.М. «Численные методы». Издательство Московского университета. 1987.

13. Самарский А.А., Гулин А.В. «Численные методы». Издательство «Наука» 1989 г.

14. Дьяконов В. «MATLAB 6: учебный курс». Издательский дом «Питер» 2001 г.

15. Paraschivoiu J. «Double Multiple Stremeamtube model with Recent Improvements» // Journal of Energy, vol.7 no.3.

16. Vittecoq P., Laneville A. «The Aerodynamic Forses for a Darrieus Rotor with Straight Blades: Wind Tunnel Measurement» // Journal of Wind Engineering and Industrial Aerodynamics vol.15 Aug-Sept. 1983 pp 381-388.

17. Самсонов В.А. Беляков Д.В. «Математическая модель движения симметричного авторотирующего тела в сопротивляющейся среде». Юбилейный сборник к 75-летию МАТИ. Изд-во МАТИ-РГТУ. 2007

18. Беляков Д.В., Самсонов В.А. «Оценка возможностей нового типа ротирующего спускающегося в воздухе объекта». XXVI Академические Чтения по Космонавтике. 2002 г.

19. Беляков Д.В. «Математическое моделирование движения ротирующего спускающегося в воздухе объекта». Труды Пятого Международного Аэрокосмического Конгресса 2006 г.

20. Беляков Д.В. «Математическая модель несимметричного авторотирующего тела в сопротивляющейся среде». ХХХШ Международная Молодежная Научная Конференция «Гагаринские Чтения» 2007 г.

21. Беляков Д.В. «Математическое моделирование движения ротирующего спускающегося в воздухе объекта». Пятый Международный Аэрокосмический Конгресс IAC 06. Посвящается 20-летию вывода в космос орбитальной станции "МИР". Полные доклады. 27-31 августа 2006 г., Москва, Россия. Электронный вид. Рег. номер.

22. Beljakov D.V. «Mathematical modelling of movement rotating object going down in air». Fifth International Aerospace Congress IAC'2006. 27-31 August, 2006, Moscow, Russia. SIP RIA, 2006. Full of lecture. Electronic form. Serial number. 2007. Л.