2. Вывод уравненений движения самолета
Уравнения движения самолета относительно инерциальной системы отсчета могут быть получены из основных теорем динамики твердого тела. Движение твердого тела описывается векторными уравнениями:
и - главный вектор и главный момент относительно центра масс количества движения твердого тела (, );
и - главный вектор и главный момент относительно центра масс внешних сил, действующих на твердое тело.
Если рассматривать самолет как твердое тело в произвольный момент времени и используя принцип “затвердевания”, то векторные уравнения количества движения и момента количества движения самолета в инерциальной системе отсчета примут вид:
- уравнение описывающее поступательное движение центра масс самолета (уравнение сил).
- уравнение описывающее вращательное или угловое движение вокруг центра масс самолета (уравнение моментов).
Где и - количество движения и момент количества движения относительно центра масс самолета как затвердевшей системы переменного состава; и - главный вектор и главный момент внешних сил, не связанных с работой двигательной установки; и - тяга двигателей и момент тяги двигателей относительно центра масс; и - главный вектор и главный момент относительно центра масс кориолисовых сил инерции.
Пренебрегая скоростью и ускорением перемещения центра масс самолета относительно его корпуса, вычисляем производную согласно принципу затвердевания:
где m - масса самолета; - абсолютная скорость его центра масс.
Пренебрегая влиянием кориолисовых и вариационных сил и моментов, связанных с движением масс топлива и газа внутри самолета, получаем:
Удобнее исследовать движение самолета, пользуясь подвижными системами координат вначале с центром масс самолета. При проектировании производной по времени от какого-либо вектора (определенного относительно системы отсчета ) на оси любой подвижной системы координат OXYZ, вращающейся с угловой скоростью относительно выбранной системы отсчета (неподвижной):
Размещено на http://www.allbest.ru/
Где - проекции вектора на оси системы ОХYZ; - их производные; - проекции угловой скорости на оси системы OXYZ.
Y
Z X
Перепишем полученную систему для вектора скорости :
Размещено на http://www.allbest.ru/
Для продольного движения самолета =
Принимая во внимание малость абсолютных величин переносной и кориолисовой сил инерции, связанных с вращением Земли:
где - вектор гравитационного ускорения; - главный вектор аэродинамических сил.
Векторное уравнение движения центра масс самолета примет вид:
где - вектор скорости движения центра масс самолета относительно Земли; - главный вектор аэродинамических сил; - сила тяжести.
Наиболее простую и удобную форму система динамических уравнений движения центра масс самолета примет, если векторное уравнение спроектировать на оси траекторной системы координат .
-кинематический угол атаки
-угол атаки ветра
- угол атаки
-траекторный угол
-угол тангажа
Применяя формулу (1) для проектирования левой части уравнения (2) и учитывая, что , получим:
для продольного движения будет отсутствовать уравнение:
Проекции аэродинамической силы на оси траекторной системы координат выражаются через проекции на скоростные оси (для продольного движения):
где -сила лобового сопротивления;- подъемная сила.
Используя матрицу направляющих косинусов между осями связанной и траекторной систем координат, проекции тяги двигателей на оси траекторной системы координат получим в следующем виде (для продольного движения):
;
Сила тяжести самолёта приложена в его центре масс, направлена по местной вертикали вниз и, следовательно, расположена в плоскости OXKYK траекторной системы координат. Её проекции на оси этой системы координат имеют вид (для продольного движения):
Тогда система динамических уравнений движения центра масс самолета в траекторной системе координат, примет вид:
Введем вектор воздушной скорости самолёта , связанный с векторами и (скорость ветра) соотношением:
Размещено на http://www.allbest.ru/
Размещено на http://www.allbest.ru/
Отсюда получим соотношение, для проекций аэродинамических сил на оси траекторной системы координат с учетом ветра (для продольного движения):
В результате система динамических уравнений движения центра масс самолета получится в виде:
Исследования движения самолёта относительно центра масс удобно выполнять, если использовать динамические уравнения в проекциях на оси связанной системы координат. При изучении углового движения самолёта так же, как и при определении траекторий центра масс, применяют в качестве системы отсчёта неинерциальную систему, связанную с Землёй.
Проектируя векторное уравнение на оси связанной системы координат и применяя формулы (1) для вычисления проекций производных по времени от вектора кинетического момента самолёта, получим систему скалярных уравнений движения самолёта относительно центра масс:
;
;
.
- проекции вектора кинетического момента самолета на связанные оси координат,
- проекции вектора абсолютной угловой скорости самолета на те же оси,
- проекции главного момента аэродинамических сил и сил тяги относительно центра масс на те же оси.
Поскольку основная плоскость ОХУ связанной системы координат является плоскостью симметрии самолета, то == 0.
Тогда система уравнений примет вид:
Для продольного движения самолета то:
Кинематические уравнения связывают между собой кинематические и геометрические характеристики поступательного движения центра масс самолёта и вращения его относительно центра масс, а также угловые скорости подвижных систем координат с параметрами движения самолёта.
Кинематическое уравнение движения центра масс самолета в векторной форме
где - радиус-вектор и вектор скорости центра масс самолета относительно рассматриваемой системы отсчета. Для получения скалярных кинематических уравнений движения центра масс найдем проекции вектора скорости центра масс самолета на оси координат, относительно которых рассматривается движение самолета. Проектируя вектор скорости на нормальные оси координат и используя таблицу направляющих косинусов, получим кинематические уравнения движения центра масс самолета:
- координата самолета в стартовых осях.
Кинематические уравнения, описывающие вращение самолета относительно нормальной системы координат, устанавливают связь между производными углов -по времени и проекциями на связанные оси вектора угловой скорости самолета относительно системы отсчета, связанной с Землей. Поскольку вращение самолета может быть представлено как изменение углов , определяющих положение самолета относительно Земли, вектор угловой скорости самолета равен геометрической сумме угловых скоростей элементарных поворотов
.
Это уравнение является кинематическим уравнением вращательного движения самолета в векторной форме. Проектируя векторы на направление связанных осей OX, OY и OZ получим:
Для продольного движения получим:
В результате получим систему уравнений продольного движения самолёта с учётом ветра:
(*)
3. Линеаризация дифференциальных уравнений движения самолёта
Линеаризация функции в окрестности значения аргументов - разложение функции в ряд Тейлора по первым членам в этой окрестности:
- опорные значения аргументов.
При линеаризации д.у. переходят от самих параметров движения к их приращениям относительно опорных значений.
Соответственно производные: .
Опорным движением является установившийся горизонтальный полёт, поэтому:
Линеаризация правых частей уравнений в системе (*):
1) в первом уравнении
2) во втором уравнении
3) в третьем уравнении
Учитывая, что ,
Выполнено следующее приближенное равенство:
так как:
Так как самолёт дозвуковой имеем:
Следующие слагаемые являются малыми и ими можно пренебречь:
Также можно пренебречь тягой по сравнению с , по сравнению с .
Получаем:
Пренебрежение слагаемыми позволяет исключить из рассмотрения уравнение для , положив .
Обозначим: , уравнения перепишутся в виде:
Из уравнений исключается и получается система уравнений:
4. Вывод передаточных функций
Достаточно рассмотреть первые два уравнения системы.
Примем, что в силу малости Wx и равны нулю.
Учитывая что , а следовательно
то проведём подстановку в первые два уравнения системы.
Проведём преобразование Лапласа, введём закон управления() и сделаем группировку относительно , и .
1) Первое уравнение:
2) Второе уравнение
Запишем получившуюся систему в матричном виде:
Найдём передаточные функции .
Найдём определители матриц в приложении MATLAB Symbolic Math и сгруппируем их относительно p:
Передаточная функция
:
5. Характеристики турбулентной атмосферы. Описание метода вычисления дисперсии
При некоторых метеорологических условиях в отдельных зонах атмосферы возникают хаотические неупорядоченные движения воздуха - турбулентность. Самолет, попадая в зону турбулентности, подвергается воздействию со стороны возмущенного потока. При этом возникает болтанка самолета - дополнительная перегрузка и угловое движение, которые при полете в спокойной атмосфере отсутствуют.
Теоретические и экспериментальные исследования привели к следующим результатам:
1. Величина пульсации скорости в пределах объема, который занимает самолет обычных размеров, существенно не меняется.
2. Пульсация скорости ветра является стационарным процессом. Компоненты этой скорости Wx, Wy, Wz являются независимыми. Статистические характеристики пульсаций скорости ветра в поперечных направлениях Wy, Wz одинаковы.
3. Спектральные плотности компонент Wx и Wy имеют следующие выражения:
При выполнении курсовой работы рекомендуется использовать частотный метод статистического анализа. Согласно этому методу дисперсия высоты определяется выражением:
где , - передаточные функции, с учетом замены . , - спектральные плотности компонент Wx и Wy.
Помимо явного способа вычисления дисперсии по формуле разработан рекуррентный алгоритм вычисления, особенно удобный для создания эффективных вычислительных программ на ЭВМ.
Согласно данному методу выражение для дисперсии можно представить в виде:
где , , А и В - полиномы с рациональными коэффициентами:
Необходимо, чтобы полином имел все нули в левой полуплоскости, а полином имел все нули в левой полуплоскости и, может быть, на мнимой оси. Кроме того, степень полинома должна быть, по крайней мере, на единицу меньше, чем степень полинома .
Если эти требования к полиномам выполнены, то дисперсия может быть вычислена по рекуррентному соотношению:
,
где с начальным условием , так же здесь ,
Параметры , , - коэффициенты полиномов и , степени которых не превосходят n.
Для получения полиномов A(p) и B(p) необходимо спектральные плотности представить в виде:
далее взять первые множители , и перемножить их с полученными передаточными функциями.
Для упрощения вычислений в выражениях сделаем замены:
Далее получим:
Чтобы использовать описанную выше методику вычисления дисперсии, необходимо произвести замену . Получим:
Используя MATLAB Symbolic Math, перемножим полученные выражения для спектральных плотностей с передаточными функциями и запишем полиномы:
Аy(p)=
Вy(p)=
6. Текст программы
6.1 Главная программа:
clear all
clc;
L=1040;%масштаб турбулентности
switch 3
case 1 %слабая турбулентность
sigmay = 0.17;
Kw=0:0.01:1;
Ka=0:0.01:1;
for i=1:size(Kw,2)
for j=1:size(Ka,2)
[A, B] = polinom_y (Kw(i), Ka(j), sigmay, L);
Disp(i, j)=dispersion(A,B,size(B,2));
end;
end;
figure(Name,Дисперсия перегрузки,NumberTitle,off)
surfc(Kw, Ka, sqrt(Disp)),grid on;
title(sigma=0.17);
xlabel(Komegaz);
ylabel(Kalfa);
zlabel(D);
colorbar
shading interp
colormap(jet)
C = min(Disp);
C1=min(C);
print -djpeg results/dispersija_turbulentnosti017.jpg
case 2 %средняя турбулентность
sigmay = 1.71;
Kw=0:0.01:1;
Ka=0:0.01:1;
for i=1:size(Kw,2)
for j=1:size(Ka,2)
[A, B] = polinom_y (Kw(i), Ka(j), sigmay, L);
Disp(i, j)=dispersion(A,B,size(B,2));
end;
end;
figure(Name,Дисперсия перегрузки,NumberTitle,off)
surfc(Kw, Ka, sqrt(Disp)),grid on;
title(sigma=1.71);
xlabel(Komegaz);
ylabel(Kalfa);
zlabel(D);
colorbar
shading interp
colormap(jet)
C = min(Disp);
C1=min(C);
print -djpeg results/dispersija_turbulentnosti171.jpg
case 3 %сильная турбулентность
sigmay = 5.99;
Kw=0:0.01:1;
Ka=0:0.01:1;
for i=1:size(Kw,2)
for j=1:size(Ka,2)
[A, B] = polinom_y (Kw(i), Ka(j), sigmay, L);
Disp(i, j)=dispersion(A,B,size(B,2));
end;
end;
figure(Name,Дисперсия перегрузки,NumberTitle,off)
surfc(Kw, Ka, sqrt(Disp)),grid on;
title(sigma=5.99);
xlabel(Komegaz);
ylabel(Kalfa);
zlabel(D);
colorbar
shading interp
colormap(jet)
C = min(Disp);
C1=min(C);
print -djpeg results/dispersija_turbulentnosti599.jpg
end
%Поиск оптимальных коэффициентов Komegaz и Kalfa
%при минимальном значении дисперсии
d_min=9999;
i_min=0;
j_min=0;
for i=1:101
for j=1:101
if(d_min>Disp(i,j))
d_min=Disp(i,j);
i_min=i;
j_min=j;
end;
end;
end;
disp(Минимальная дисперсия Disp= )
d_min
disp(Komegaz= )
Kw(i_min)
disp(Kalfa= )
Ka(j_min)
6.2 Подпрограмма расчета полиномов Ay(p) и By(p):
function [A, B] = polinom_y(Kw, Ka, sigmay, L)
%%%геометрические и физические характеристики%%%
S=160;%площадь крыла
ba=4.7;%САХ
m = 86000;%масса самолёта
Iz=6100000;%момент инерции
xt=0.4;%центровка самолёта
Lkab=15.8;%плечо кабины экипажа
g=9.81;
%%%станартная атмосфера%%%
M = 0.695;%наивыгоднейшее число M самолёта
a=303.9;%скорость звука
ro=0.467;%плотность воздуха
Vnv=M*a;%наивыгоднейшая скорость самолёта
V0=Vnv;
q=(ro*Vnv^2)/2;%скоростной напор
%%%аэродинамические характеристики%%%
mzwz=-10.4;
mzd=-1.32;
Cyalfa=5.61;
xf=0.78;%фокус самолёта
mza=(xt-xf)*Cyalfa;
Cy_gp=m*g/(q*S);
alfa0=Cy_gp/Cyalfa;
Kz=q*S*ba/Iz;
Mzwz=mzwz*Kz*(ba/V0);
Mzd=mzd*Kz;
Mza=mza*Kz;
d1=(1)/(g*V0);
a2=alfa0*V0;
a1=-((-g+Mzwz*V0*alfa0+Kw*Mzd*V0*alfa0)/(1));
a0=-((Mzwz*g+Kw*Mzd*g+Mza*V0*alfa0+Ka*Mzd*V0*alfa0)/(1));
b2=(V0*g+Mza*V0*alfa0*Lkab)*d1;
b1=-(Mzwz*V0*g+Ka*Mzd*g*Lkab+Kw*Mzd*V0*g)*d1;
b0=-Ka*Mzd;
k1=((sigmay^2)*L)/(2*pi*V0);
k2=L/V0;
%матрица коэффициентов полинома В(р)
B(1)=3^(1/2)*b2*k1^(1/2)*k2;
B(2)=((k1)^(1/2))*(b2+3^(1/2)*b1*k2);
B(3)=((k1)^(1/2))*(b1+3^(1/2)*b0*k2);
B(4)=k1^(1/2)*(b0);
%матрица коэффициентов полинома A(р)
A(1)=a2*k2^2;
A(2)=a1*(k2^2)+2*a2*k2;
A(3)=a0*(k2^2)+2*a1*k2+a2;
A(4)=a1+2*a0*k2;
A(5)=a0;
6.3 Подпрограмма расчета дисперсии:
function res = dispersion(A, B, N)
pi = 3.141592654;
a = A; b = B;
c = 0;
ier = 0;
if a(1) <= 0
res = NaN;
ier = 1;
return;
end;
for k = 1:1:N
if a(k + 1) > 0
alf = a(k) / a(k + 1);
bet = b(k) / a(k + 1);
c = c + bet^2 / alf;
k1 = k + 2;
if (k1-N) <= 0
for i = k1:2:N;
a(i) = a(i) - alf * a(i + 1);
b(i) = b(i) - bet * a(i + 1);
end;
end;
else
res = NaN;
ier = 1;
return;
end;
end;
res = pi * c ;
7. Результаты вычислений и выводы
График зависимости дисперсии перегрузки от коэффициентов усиления системы управления Komegaz и Kalfa при sigma=0.17:
График зависимости дисперсии перегрузки от коэффициентов усиления системы управления Komegaz и Kalfa при sigma=1.71:
График зависимости дисперсии перегрузки от коэффициентов усиления системы управления Komegaz и Kalfa при sigma=5.99:
Из графиков видно, что при изменении коэффициентов усиления системы управления Komegaz и Kalfa дисперсия перегрузки изменяется и имеет ярко выраженный минимум. Из этого следует, что введение данного закона управления при полёте в турбулентной атмосфере оказывает положительное влияние.
Оптимальные значения коэффициентов усиления:
Komegaz=
ans =
0.1600
Kalf=
ans =
0.2800
8. Список используемой литературы
самолет атмосферная турбулентность
1. Овчаренко В.Н., Павлов К.А. Методические указания к курсовой работе по теме «Статистическая динамика», М.:Изд-во МАИ,1993.
2. «Аэромеханика самолета» Бочкарев А.Ф. , М.:Машиностроение 1985г.
3. «Динамика полета в неспокойной атмосфере» Доброленский Ю.И, М.:Машиностроение 1969г.
4. «Управление полётом самолётов» Гуськов Ю.П. , Загайнов Г.И, М.:Машиностроение 1991г.
5. «MATLAB 7.*/R2006/R2007: Самоучитель» Дьяконов В. П., М.: «ДМК-Пресс» 2008г.