CC BY
27
International Journal of Open Information Technologies ISSN: 2307-8162 vol. 3, no. 2, 2015
Оптимизация спуска осесимметричного авторотирующего тела
в воздушной среде
Д.В. Беляков
Аннотация - Рассматриваются вопросы моделирования спуска авторотирующих тел сложной конфигурации в квазистатической среде. Показано существование режима авторотации с минимальной вертикальной составляющей скорости, что позволяет применять конструкцию как систему спуска. Проводится оптимизации угла отклонения скорости центра масс авторотирующего тела от вертикали. Предложен численный алгоритм оптимизации с управлением в виде гармонического закона.
Ключевые слова - тело сложной конфигурации, стационарный режим, осреднение, закон управления, оптимизация
I. Введение
В задачах динамики тела, взаимодействующего со средой, существует проблема моделирования спуска быстро вращающихся тел разной конфигурации с минимальной вертикальной составляющей скорости. Эта проблема является актуальной как в математическом, так и в прикладном направлениях. Построенные модели обычно используют в качестве систем спуска, или торможения при движении в атмосфере. Необходимость решения таких задач существует и в авиаракетостроении. Также актуальное значение имеют задачи оптимизации и управления движением рассматриваемых тел.
II. Построение математической модели рассматриваемого тела
Рассмотрим движение в воздушной среде тела сложной конфигурации, состоящего из стержня и двух
параллельных пластинок. Плоскости пластинок образуют угол S с плоскостью, ортогональной стержню (Рис. 1). Будем считать, что движение происходит в вертикальной плоскости.
При создании модели воздействия среды на тело используется гипотеза о квазистационарном обтекании пластинок средой [3], [4], [5], [6], [12] Согласно этой гипотезе сила воздействия среды на каждую пластинку характеризуется скоростью некоторой ее точки, которая называется центром давления. Помимо упомянутой используются вихревая [7] и многотрубчатая [9] модели взаимодействия пластинок со средой. Активно исследуется влияние неоднородности потока на
Статья получена 10 января 2015.
Д. В. Беляков, к.т.н., доцент, МАТИ-РГТУ им. К.Э. Циолковского, г. Москва (e mail: [email protected])
аэродинамические нагрузки [8] и эффект срыва потока. В рассматриваемой модели предполагается, что поперечные размеры пластинок намного меньше длины стержня. При выполнении этого условия, центры давления пластинок точки А и В можно считать неподвижными относительно пластинок. Считается также, что среда не оказывает никакого влияния на стержень. Аэродинамические силы, действующие на каждую пластинку, разложим на две составляющие: силы сопротивления SA, SB,
направленные против абсолютных скоростей у ,у центров давления, и подъемные силы PA, PB,
направленные им ортогонально.
При этом величины аэродинамических сил равны:
|SA |= s(a+S)Vl = 0.5pacx(a+S)V2A,
| PA |= p(a+ S)V2A = 0.5pacy (a+S)V2 | SB |= s(fi+S)V2 = 0.5pacx (P+S)V2B,
| PB |= p(P+S)V2B = 0.5pacy (P+S)V2B
Где а, в - углы атаки между векторами окружных скоростей (VA0, VBO) точек A, B и векторами у, у, p, S - аэродинамические функции углов атаки, Cx , Су -безразмерные аэродинамические функции, р - плотность воздуха, (Г - площадь одной пластинки. Зависимости cx (а) и cy (а) определены из продувок прямоугольных
пластинок с заданным удлинением в аэродинамической трубе и являются экспериментальными данными [6]. Их
11
International Journal of Open Information Technologies ISSN: 2307-8162 vol. 3, no. 2, 2015
типичный вид (для удлинения X — 8) представлен на рисунке 2.
В качестве обобщенных координат, определяющих положение тела, введем координаты x , y центра О масс, совпадающего с серединой стержня АВ, и угол / отклонения стержня АВ от вертикали. Для описания распределения скоростей точек нашего тела, зададим
величину вектора абсолютной скорости центра масс V ,
угол Y отклонения вектора V от вертикали, угол в отклонения стержня AB от вектора абсолютной скорости центра масс V и абсолютную угловую скорость стержня
a.
а
Рис. 2 Аэродинамические функции: а) подъемной силы; б )силы сопротивления
В системе координат, связанной с вектором V (см. рис. 1), для векторов VAO ,VA, VB имеем представление: VAO — {-/bsin в, racos в}, VA — {V - rasin в, racos в}
VB — {V + rasin в, - racos в}
Следовательно, силы сопротивления можно определить как:
SA — s(a+S)VA{rюsinв-V, -rоcosв},
SB — s(fi + S)VB {-V - rasin в, racos в} подъемные силы будут иметь вид:
PA — p(a+S)VA{racosв, rasi^-V},
PB — р(в+S)VB{-racosв, -V-rasinв}
Вычислим момент аэродинамических сил, приложенных к пластинкам относительно центра масс.
М0 — r(VA2(p(a+ S) sin a- s (a+S) cos a) + V'2 (р(в + S)sin в~s(P + S) cosp)ez
Для того, чтобы составить уравнения движения тела, воспользуемся теоремой о движении центра масс в проекциях на направление вектора V и ему
ортогональное, теоремой об изменении кинетического момента и теоремой о сложении скоростей. Получим:
mV—p(a+S)VAracosв- p(в+S)VBracosв+s(a+S)VA(rоsinв-V)--s(fi+SVB(ramв+V) +mgcos у
nV'У—-s{a+S)^Arъasв+s(в+SVBrъasв+p(a+SVA(ramв-V) - (1)
-p(в+S)VB(rшnв+V)-mgsiny
Но— r(V2(p(a+S)sina-s (a+S)cosa)+V;(p(в+Ssinв-s(в+Scosв) в + у—a
Соотношения, определяющие величины va ,Vs ,a, в через псевдоскорости V, в, a, имеют вид:
VA sina — -Vcosв,VA cosa— ra-Vsinв (2)
VB sin в — V cos в, VB cos в — ra+V sin в
После того, как мы интегрирования системы уравнений
(1) -(2), мы можем окончательно определить положение тела при помощи интегрирования кинематических соотношений (3)
x — V cos у, y — V sin у и вычисления / — Y + в . (3)
Таким образом, математическая модель определяется замкнутой системой уравнений
(1)-(3) , которая полностью описывает движение
рассматриваемого тела
III. ПОСТАНОВКА ЗАДАЧИ УПРАВЛЕНИЯ
В [1], [2] при помощи метода осреднения были получены формулы стационарного режима авторотации. Математические основы метода осреднения описаны в [10]. Пластинки устанавливались на постоянный малый установочный угол S относительно нормали к стержню. Абсолютная величина угла S ограничена некоторым значением S*. Установившийся угол спуска Y0 на режиме авторотации пропорционален величине угла перекоса пластинки S. Например, в случае прямоугольных пластинок с удлинением X — 8, и установочном угле S — 50, абсолютная величина
соответствующего ему стационарного значения
I Y01—13.86760, а численное решение дает еще меньшее значение. Возникает вопрос: каким образом можно
увеличить Y0 за счет изменения угла S? Аналогичная задача рассматривалась в [12]. В настоящей статье будет рассмотрена задача управления установочным углом S для увеличения стационарного значения 0 , т.е. нужно выбрать закон изменения угла S такой, чтобы достигался максимум стационарного значения 0 . В этом случае
установочный угол зависит от $ , а уравнения движения тела (1)-(3) не изменятся. При выборе управляющего воздействия S($) будем считать, что пластинки
совершают гармонические колебания с периодом 2п около некоторого положения S0, так что
S($) — S0 + B sin$+ C cos$, и задача сводится к выбору наилучших значений B, C, S0 .
12
International Journal of Open Information Technologies ISSN: 2307-8162 vol. 3, no. 2, 2015
IV. ОСРЕДНЕНИЕ УРАВНЕНИЙ ДВИЖЕНИЯ Будем считать, что тело совершает спуск в режиме авторотации и закон изменения угла 8 одинаков для обеих пластинок. Уравнения движения (1) с учетом безразмерных аэродинамических функций cx , Cy будут иметь вид:
mV =1po[cy(a+8)VAr^cy&e-cy(e+8)Vslr^cy&e+Cx(a+8)VA(rmsaie-V)-Cx(e+8)VB(rmsaie+V)]+mgco&Y mVy= 1ptT[-Cx(a+8)V1rfficos^+Cx(e+8)VBrfficos^+Cy(a+8)V1(rfflsin^-V)-Cy(e+8)VB(rfflsin^+y)]-mgsinY Ja=1por[VA2(Cy(a+8)sma-Cx (a+8)ao.d)+Vl(Cy(fi+8)&ne-Cx(e+8)ccy&pj] в+у=a
Аэродинамические функции Cx, Cy разложим в ряд учитывая, что а, в ^ 0
Cx (а+8) = Cx8 + Сх8а, Cy (а + 8) = Cy8 + Су8а,
Cx (в + 8) = Cx8 + ХР’ Cy (в+8) = Cy8 + ХР
где: Cx8 = Cx 8 X = CX 8 Cy8 = Cy 8 Cy8 = C'y (8)
Левые части кинематических соотношений (2) линеаризуем при а,в ^ 0
VAa = -V cos в, VBfi = V cos в,VA = ra-V sin в, VB = ra + V sin£ Пользуясь этими соотношениями, сделаем в уравнениях (1) , переход от переменных a,p,VAVB к переменным
e,V ,ю и получим систему уравнений (4):
<
J a = p<Jr[-Cxj(ra>)2 + V2 (Cy8 cos2 в- Cx8 sin2 в) + -V (Cy8 - Cx8) sin 2в] mV = -p&rVa[C'y8 cos2 в + Cx8(1 + sin2 в) + ~(Cy8 + cX8)sin2в] + mg cos у mV у = parVa[(a'78 cos2 в- Cy8(1 + sin2 в) + - cX8)sin2в] - mg sin у
(4)
“Коэффициенты C s,C'g,C C'g связаны с управляющим воздействием 8 следующим образом:
Cx8 = Cx0 + Cx28 ’ C8 = 2Cx28, Cy8 = C'y08, C'y8 = C'y0 , где Cx0 = Cx (0). Cx2 = ^2C’(0)' Cy0 = C'y (0)
Подставим эти соотношения в уравнения (4) и проведем осреднение их правых частей по в. Осредненная система имеет вид (5):
'Ja = p<Tr[-(ra)2(Cx0 + Cx2(r02 +1B2 +1C2)) + V2(1(C,y0 -Cx„) +
2 2
3 ,
1,
<
+Cx2(- Y2 +- B2 + - C2)]
2
8 8
1 3 3,7,5,
mV = -porVa-Cy0 ^-Cx0 + Cx2(- Yc +—B +~C )] + mgcosy
2 2 2 8 8
3 , 3 1
„mVy = parVO[(- 2 Cy0 + c^^ + 44 BC)] - mg smy
(5)
Система (5) имеет стационарное решение (6):
a= в + у
1
2mgCosy) 3 1 Cy0 - Cx0 - Cx2(8j2 + 4B2 + 4 C2)
pa(Cy0 + 3Cx0 + 3CxA2 + Cx2(^7 В2 + J4 C2)) ^ 3 1 2c + Cx2(8j2 + - B2 +1C2) x0 x2K 0 4 4
2mgCosy0 3 1 2Cx0 + Cx2(8J2 + 4 -®2 + 4 c2)
po(Cy0 + 3Cx0 + 3Cx2802 + Cx2 (4 B2 + - C2))^
'y0 Vx0 '■'xl
2Cy2(8„ + 4BC) - 3Cy„8„ -----------------------7------5—]
C7 + 3Cx0 + Cx2(38,2 + - B2 + - C2)
Y = arCtg[-
(6)
r
Из формул (6) при B = 0 и C = 0 получаются формулы режима авторотации. Задача оптимизации состоит в том,
чтобы найти такие коэффициенты B, C и значение установочного угла 8 = 80 чтобы функция
у= у( B, C,8) достигала максимального значения.
V. ОПТИМИЗАЦИЯ Найдем стационарные точки функции (7) :
z(x, y) =
80 (2Cx2 - ЗсУр) + 2 Cx2xy = a + hxy
Cy 0 + 3Cx0 + ЗГ2802 + 7 x2 + 5 y2 C d X + ^ y
(7)
.Для определения стационарных точек, будем иметь систему уравнений:
b1 y(a1 + d1 xV + l1 y2) = 2d1 x[a1 + b1 xy] (8)
b x(a1 + d1 x +11 y2) = 2l1 y[a1 + b1 xy]
Следовательно, x = 0 и y = 0 являются ее решениями. Предположим, что x Ф 0 и разделим первое уравнение на
y d 1 x
вт о р о е . Бу д е м и м е т ь : — =-, от куда:
x l y
y =±
d1 /7
—x = ±k1 x. Подставляя y =J—x и значения
l1 V5
коэффициентов a1, b1, C1, d1, l1 в первое уравнение системы (4.2), получим:
Cx2^/5(CyC + 3Cx0 + 3Cx2 80 + 2 x ) = 7[80 (2Cx2 - 3Cy0 ) + 2 Cx2x ]
, откуда следует, что у этого уравнения отсутствуют
7
решения. Подставляя y = -у—x в первое уравнение системы (4.1), получим:
Cx2^4(cyc + 3Cx0 + 3Cx2 80 + 2 x ^ = 7[80(2Cx2 - 3Cy0 ) - 2^5"Cx2x ]
Окончательно получим:
80(14Cx2 - 21<0) - Cx^J'-(CУ
x2 =
y0' '-'x2\l - VL'y0 + 3Cx0 + 3Cx280 )
< 0
XC-
Таким образом, единственная стационарная точка
исследуемой функции, существует только в точке {0,0}. Найдем частные производные второго порядка в найденной стационарной точке:
13
International Journal of Open Information Technologies ISSN: 2307-8162 vol. 3, no. 2, 2015
„ 2d, a
Zxx (0,0) = 11
7(3<0 - 2cx2)8
c1 2(c'y 0 + 3cx0 + cx2 8 )2
Zyy (0,0) = -
5(3cy0 - 2cx2)£
2l1a1
c1 2(cy0 + 3cx0 + cx238 )2
b
Zxy (0,0) = ^ =
c1 4(cy0 + 3cx0 + cx2 38 )
Составим дискриминант:
D = Zxx (0,0)Zyy (0,0) -(Zxy (0,0))2 =-
4d1l1 a1 b1
,2
c
c
c
x2
При выполнении условия D > 0 в точке {0,0} функция (7) будет иметь минимум, если D < 0, то он превращается в седловую точку. Найдем значения параметра 8 = д0, при которых выполняется неравенство D > 0 .
4d1l1a12 b21
----2- > 0 ^ 4d1l1 af - b12 cf > 0 ^ 4d1l1a12 > b12 cf ^
^35((3cy0 - 2cx2)282 > (cy0 + 3cx0 + 3cx2d2)
y0
4
c
c
После очевидных преобразований, получим следующее неравенство:
то максимум z(x) равен a1 = (2cx2 - 3c^0 )8. Если
a1 < 0 , или 0 < 8 < 81, то максимум достигается при
x .
Г Ьхк
lim z( x)
x^<x
2
(d1 + l1k2) к=
, = b1(d1 + lxk ) - 2l1b1 Г = 0 k = (d1 +11 к2)2 = ’
V
di 4l Ф=^ - 9.680
b1k
(d1 + l1k2) ’ откуда
Рассмотрим теперь случай минимума. В этом случае надо найти максимальное значение функции (7) на
окружности + Б2 =8* - д0 = R. Сделаем замену
переменных: A = R cos <р, Б = R sin <р, подставим ее в функцию (7), получим:
z(q>) =
a1 + — b1 R2 sin 2^
c1 + d1 R2 cos2 ^ +11R2 sin2 ^
, где z(^) является
нечетной 2п периодической функцией. Достаточно найти ее максимум на отрезке [0 п]. Изобразим график z(^) для прямоугольных пластинок с удлинением Л = 8 (рис. 3). Из рисунка 3 убеждаемся в том, что абсолютное значение максимума z(0) близко к нулю.
3cx2282 - V35(3<0 - 2cx2)1 8 11 +cx2 (S0 + 3cx0 ) < 0 (9)
Правая часть (9) обращается в нуль в следующих точках:
§ _л/35(3с;0 -2сх2)-ф5(3с'у0-2сх2)2 -12сх2(с'у0 +3сх0] ^ ^
6cr
82 =-
л/35 (3cy0 - 2cx2 ) +yj35(3c'y0 - 2cx2 )2 - 12cx2(cy0 + 3cx0 )
"6с
X
В случае прямоугольных пластинок с удлинением Л = 8: 8 -1.50,82 -9.530
О 044 I-----------1-----------1-----------1-----------1------------1-----------1-----------1-----------1-----------1-----------1
0 042 --- ......\..........-
0 04 — ■ ■ ■ ■ I ■ I ■ ■ ■ .......г..........-
0 035 --- : : : \ / \ \ \ ; ; ......I..........-
0 034 --- ! ! : : / : \ : : : ......i..........-
0 032 -..........i...........i............:.......................j.......................-i...........i-..........i..........-
0 03 I-----------1----------1------------1-----------1------------1-----------1-----------1-----------1-----------1-----------
-10 -8 -5 -4 -2 0 2 4 5 8 10
Рис.3 график функции z(0)
Решив неравенство (9), убеждаемся в том, что в стационарной точке {0, 0} возможны следующие случаи:
1. Если |8 |< 81 , то D < 0 и функция (7) будет иметь седловую точку.
2. Если 8 > 81 , то D > 0, A > 0 и (7) будет иметь минимум.
3. При 8 <- 81, D > 0, A < 0 и (7) будет иметь максимум.
Рассмотрим сначала случай точки типа «седло». В этом случае существует прямая y = kx, на которой
Мы видим, что, при спуске авторотирующего тела, оптимизация средней боковой силы по первому приближению с управляющим воздействием 8(8) в
виде гармонического закона не приводит к существенному увеличению абсолютного значения угла отклонения от вертикали, а рассмотрение высших приближений связано с значительным объемом вычислений. Построение закона управления иным способом проводилось в [11].
Таким образом, аналитический поиск решений вызывает затруднения. Значение максимума и коэффициентов Б, C , можно определить численно.
достигается максимум.
z( x) =
a1 + b1 kx2
c1 + ( d1 + l1 k ) x
Функция (7) примет вид: Если a1 > 0, или -81 <8 < 0 ,
VI. ОПТИМИЗАЦИЯ УГЛА ОТКЛОНЕНИЯ ОТ ВЕРТИКАЛИ ПРИ ПОМОЩИ ЧИСЛЕННЫХ ИССЛЕДОВАНИЙ
14
International Journal of Open Information Technologies ISSN: 2307-8162 vol. 3, no. 2, 2015
Задача поиска максимума функции у=у0 (B, C) может
быть сведена к поиску минимума функции f =-у (B, C). Для минимизации f будем использовать встроенные
функции Matlab fminu и fminsearth. Эти функции реализуют метод скорейшего градиентного спуска и симплекс-метод Нелдера - Нида.
Рассмотрим решение такой задачи: найти минимум функции
f = -arctg[
2Cx2 ( Y0 + 4 BC) - 3c'y0
7 5
СУ0 + 3cx0 + cX2(3Y + 4 B2 + 4 C2)
Введем значения аэродинамических функций и угла перекоса 8 = 50. Такому установочному углу соответствует стационарное значение | у |= 13.8676°- В
любом вычислительном процессе, использующем
итерационные методы, необходимо указать начальную точку. Пусть в нашем случае этой точкой будет (0Л;0Л). Теперь необходимо вызвать функцию fminu .
del=5
del=del*(pi/180)
hh=(pi*5)/180
s=[1.4 1.38 1.35 1.3 1.25 1.17 1.09 0.99 0.89 0.79 0.71 0.61 0.5 0.4 0.28 0.2 0.11 0.02 0.01 0.02 0.11 0.2 0.28 0.4 0.5 0.61 0.71 0.79 0.89 0.99 1.09 1.17 1.25 1.3 1.35 1.38 1.4]; p=[0 -0.17 -0.28 -0.4 -0.5 -0.6 -0.65 -0.72 -0.77 -0.81 -0.87 -0.9 -0.87 -0.79 -0.73 -0.75 -0.78 -0.4 0 0.4 0.78 0.75 0.73 0.79 0.87 0.9 0.87 0.81 0.77 0.72 0.65 0.6 0.5 0.4 0.28 0.17 0]; d1CY0=p(:,20)-p(:,19))/hh CX2=0.5*(s(:,20)-s(:,19))/(hhA2)
CX0=s(:,19)
CXDEL=CX0+CX2*delA2 CYDEL=d1 CY0*del d1CXDEL=CX2*del
ga=(d1CXDEL-3*CYDEL)/(d1CY0+3*CXDEL)
%ga=180*ga/pi
A=atan(ga)
global A d1CY0 CX2 CX0 d1CXDEL x y p s x=[2 4]
xminvector = fminu(@minfun1,x);
после завершения работы программы мы видим, что B=-0.04669 C=-0.02446 и | y0\= 0.6366 рад. или I y0 = 36.4745 1
VII. ЗАКЛЮЧЕНИЕ
Таким образом, мы ввели управляющее воздействие на установочный угол пластинок ^ = S0 + B sin^+ C costf,
провели осреднение правых частей уравнений движения, с помощью численных методов провели оптимизацию угла отклонения от вертикали при движении в режиме авторотации с управляющим воздействием S($), что позволило увеличить угол отклонения в 3 раза.
сопротивляющейся среде». Издательство «МЭИ», журнал «Вестник МЭИ», выпуск № 4 2006 г., стр. 5-10.
[2] Беляков Д.В. "Исследование и особенности
математической модели движения несимметричного авторотирующего тела в квазистатической среде". Издательство «Новые технологии», журнал " Мехатроника, Автоматизация, Управление". Выпуск № 11. 2007 г., стр. 20-24
[3] Локшин Б.Я. , Привалов В.А., Самсонов В.А. “Введение в задачу о движении точки и тела в сопротивляющейся среде”. Издательство Московского университета. 1992.
[4] Паршин Д.Е. Самсонов В.А. «Качественный анализ в задаче о движении аэродинамического маятника.»1992г. (МГУ НИИмеханики, отчет 419).
[5] Б.Я. Локшин, В.А. Привалов, В.А. Самсонов «Качественный анализ движения твердого тела в сопротивляющейся среде».Шестой Всесоюзный съезд по теоретической и прикладной механике. Ташкент, 1986 г.
[6] Табачников В.Г. «Стационарные характеристики
крыльев на малых скоростях во всем диапазоне углов атаки.» Труды ЦАГИ 1974 г. выпуск 1621
[7] Strickland J.N. , Smith T. and Sun K. ,”A vortex Model of the Darrious Turbine: An analitycal and experimental stady”, Final Report submitted to Sandia Laboratories on contract # 13-5602, 1982
[8] Г.И. Шоломович, В.П. Карликов, О.В. Зверев «Исследование по оптимизации ортогональных энергоагрегатов для использования энергии воздушных и водных течений. Оценка влияния кривизны и неоднородности потока на на аэродинамические характеристики профиля»,-М.,1990 г.-33 л.(МГУ. НИИмеханики, отчет 3910)
[9] Paraschivoiu J., “Double Multiple Stremeamtube model with Recent Improvements”, Journal of Energy, vol.7 no.3
[10] В.Ф. Журавлев Д.М. Климов «Прикладные методы в теории колебаний»
Издательство «Наука» 1988 г.
[11] Самсонов В.А. Беляков Д.В. «Оптимизация
движения спускающегося авторотируюшего тела». Седьмой Международный Аэрокосмический Конгресс IAC'07. Полные доклады. 27-31 августа 2012 г., Москва, Россия. Электронный вид. Регистрационный номер: ISBN 7-85312-056-7
[12] Имитационное моделирование спуска
осесимметричного авторотирующего тела в квазистатической среде. Девятая Международная научнопрактическая конференция «Современные
информационные технологии и ИТ образование» 13-15 ноября 2014 г. Сборник трудов стр. 856-864
Библиография
[1] Беляков Д.В., Самсонов В.А., Филиппов В.В. «Исследование движения несимметричного тела в
15
Optimization of descent axisymmetric windmilling body in the air
D.V. Belyakov
Abstract - They Are Considered questions of modeling of the lowering autorotations bodyes complex deskside in quasistatic to ambience. Existence of the mode autorotations is shown with minimum vertical forming velocities that allows to use the design as system of the lowering. It Is Conducted optimization of the corner of the deflection to velocities of the centre of the masses autorotations bodies from vertical. It Is Offered numerical algorithm to optimization with management in the manner of harmonic law.
Keywords - the body is a complex configuration, the stationary mode, averaging, law management, optimization
16
