98946

Анализ влияния атмосферной турбулентности на продольное движение самолета (Прототип Ягуар, фронтовой истребитель)

Отчет о прохождении практики

Астрономия и авиация

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

Русский

2016-07-17

1011 KB

0 чел.

Отчёт

по курсовой работе по дисциплине

“Статистическая динамика”

Тема: «Анализ влияния атмосферной турбулентности на продольное движение самолета»

(Прототип “Ягуар”, фронтовой истребитель)

Содержание.

  1. Аннотация………………………………………………………………………………………..3
  2. Задание…………………………………………………………………………………………...4
  3. Исходные данные к работе……………………………………………………………………...5
  4. Вывод уравнений, описывающих продольное движение самолета при наличии ветра……6
  5. Линеаризация  уравнений  движения . ………………………………………………………..11
  6. Вывод передаточных функций……………………………………………………………….. 14
  7. Характеристики турбулентной атмосферы .Метод вычисления дисперсии………………..17
  8. Программная реализация………………………………………………………………………22
  9. Результаты расчетов……………………………………………………………………………35
  10. Список используемой литературы…………………………………………………………… 39

Аннотация.

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

При этом решаются две основные задачи:

  1. Расчет влияния атмосферной турбулентности на характеристики движения самолета

     без выключенной САУ;

  1. Расчет влияния САУ на характеристики движения в условиях атмосферной турбулентности.

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

Задание на курсовую работу.

1. Разработать математическую модель полёта самолёта в турбулентной атмосфере.

а) Вывести линеаризованные уравнения продольного движения самолета при наличии ветра.

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

в) Из системы ДУ получить передаточную функцию, описывающую изменение выбранной переменной в зависимости  от  скорости ветра.

2. Вычислить зависимости дисперсии выходного параметра от режима полёта, параметров

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

3. Провести анализ полученных результатов.

Исходные данные

S = 25 м2– площадь несущей поверхности.

ba = 2.8 м– средняя аэродинамическая хорда.

Iz = 2.05 * 105 кг * м2– момент инерции относительно осиZ.

– относительное положение центра тяжести.

m0 = 15000 кг – взлетная масса самолета.

g = 9.81 м / с2 – ускорение свободного падения.

=6 м – расстояние от центра тяжести до кабины

P – атмосферное давление.

a – скорость звука.

– плотность воздуха.

H, м

2000

4000

6000

8000

10000

12000

P(Н), Па

79490

61660

47210

35650

26290

19390

a(Н), м /c

332.7

324.7

316.6

308.2

299.6

295.2

(Н), кг / м3

1.01

0.819

0.66

0.526

0.414

0.312

– производная коэффициента подъемной силы по углу атаки.

– относительное положение фокуса.

– производная коэффициента момента тангажа по отклонению руля высоты.

– производная коэффициента момента тангажа по угловой скорости.

– производная коэффициента момента тангажа по производной угла атаки.

М

0.5

0.6

0.7

0.8

0.9

1

1.1

1.2

1.3

, 1/рад

3.6

4

4

4.1

4.2

4.1

4.5

4.7

4.8

0.36

0.4

0.43

0.425

0.42

0.45

0.5

0.56

0.6

-1.126

-1.247

1.249

-1.248

-1.247

-1.248

-1.251

-1.4

-1.35

-3.3

-3.3

-3.3

-3.3

-3.3

-3.4

-3.6

-4

-3.75

-1.06

-1.04

-1.02

-0.98

-0.92

-0.87

-0.87

-0.62

-0.55

L = 1000 м – масштаб турбулентности.

= 5 – дисперсия горизонтальной составляющей порывов ветра.

= 5 – дисперсия вертикальной составляющей порывов ветра.

коэффициент усиления по перегрузке.

коэффициент усиления по угловой скорости.

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

0.5

1

1.5

2

2.5

3

3.5

4

4.5

5

Вывод уравнений движения самолета

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

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

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

;               ;

где  и - главный вектор и главный момент относительно центра масс количества движения твердого тела (, );

и   -  главный вектор и главный момент относительно центра масс внешних сил, действующих на твердое тело.

Если рассматривать самолет как твердое тело в произвольный момент времени, то к нему будут приложены:

  • внешние силы, действующие на систему;
  • сила тяги двигателя;
  • внутренние кориолисовые силы инерции

Тогда уравнения движения самолета примут вид:

;                      ;

Удобнее исследовать движение самолета, пользуясь подвижными системами координат в начале с центром масс самолета. При проектировании производной по времени от какого-либо вектора  на оси любой подвижной системы координатOxyz, вращающейся с угловой скоростью относительно выбранной системы отсчета (неподвижной), должны быть применены известные из векторного анализа формулы:

    

1) Движение центра масс самолета.

Пренебрегая скоростью и ускорением перемещения центра масс самолета относительно его корпуса, производная от количества движения по времени будет равна:

Учитывая правила векторного произведения  (*) и разделяя полученные уравнения на проекции по осямX,Y,Z, получим систему динамических уравнений движения в проекциях на оси системы координат, помещенных в центр масс самолета:

      (1)

Векторное уравнение движения центра масс с учетом того, что

- главный вектор аэродинамических  сил, приложенный в центре масс самолета;

- сила тяжести;

примет вид:

(2)

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

Применяя формулу (1) для проектирования левой части уравнения (2) и учитывая, что ; , получим:

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

Подставляем полученные проекции в первые два уравнения системы (3):

Допущения:

- нас интересует продольное движение самолета ()

- считаем самолет жестким телом, не учитываем кориолисовые силы инерции

- не учитываем слагаемые, связанные с вращением Земли;

- считаем массу самолета постоянной.

Тогда динамические уравнения движения центра масс в траекторной системе:

2) Уравнения, описывающее вращательное или угловое движение вокруг центра масс.

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

- проекции вектора кинетического момента самолета на связанные оси координат;

-  проекции вектора угловой скорости самолета относительно Земли на те же оси;

- проекции главного момента аэродинамических сил и сил тяги относительно центра масс на те же оси.

Проекции вектора кинетического момента

Тогда система уравнений примет вид (без учета угловой скорости суточного вращения Земли, угловой скорости самолета относительно нормальной системы координат и угловой скорости, возникающей из-за кривизны поверхности ):

Допущения:

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

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

В итоге получим:

3) Кинематические уравнения.

Кинематическое уравнение движения центра масс самолета в векторной форме:

,

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

- координата самолета в стартовых осях;  Н-высота полета.

В поставленной задаче используется только второе уравнение:

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

Это уравнение является кинематическим уравнением вращательного движения самолета в векторной форме. Проектируя векторы  на направление связанных осейOX, OY иOZ получим

Допущения:

- нас интересует продольное движение

4) Система дифференциальных уравнений, описывающих продольное движение самолета:

5) Модель продольного возмущенного движения самолета при наличии ветровых воздействий.

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

- угол тангажа и скоростной угол тангажа

- вектор воздушной скорости

- вектор скорости ветра в турбулентной атмосфере

Wx,Wy  - проекция ветра на осиOXa,OYa

αK - кинематический угол атаки

Поскольку аэродинамические коэффициенты обычно определяются в скоростной СК, то удобно спроектировать аэродинамические силы на нее:

Система дифференциальных уравнений, описывающих продольное возмущенное движение самолета:

Линеаризация математической модели

Линеаризация уравнений производится на основе принципа малых возмущений и с использованием разложения в ряд Тейлора нелинейных составляющих сил и моментов при сохранении величин не более 1-го порядка малости. В качестве опорного движения принимается  горизонтальный полет с постоянной скоростью в спокойной атмосфере. При этом в опорном движении

В опорном движении  имеет место следующая система:

Не учитываем угол установки и изменение тяги двигателя, пренебрегаем механизацией самолета ().

В соответствии с принципом малых возмущений:

Имея в виду, что

проведём линеаризацию правых частей уравнений в системе (**):

1) в первом уравнении

2) во втором уравнении

3) в третьем уравнении

4) в четвертом уравнении

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

где .

Следует иметь в виду:

.

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

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

Введя переменную  вместо переменной  и положив , получим:

Систему уравнений возмущенного движения  можно переписать:

где .

Вывод передаточной функции

Выходом является

Достаточно рассмотреть первые два уравнения.

Перегрузка в кабине находится из уравнений:

В итоге получим окончательное уравнение для нахождения

                                                                                                                       (1)

   (2)

Подставляем в (1):

Анологично (2) приводятся к переменным ,

Подставляя в (2):

Введем оператор Лапласа  и группируем коэффициенты при  и :

Перепишем систему в матричной форме:

Найдем определители матриц и сгруппируем результат относительноp в программеMathCad:

=

1=

2=

Пусть

;

;

;

;

;

;

;

.

Тогда:

Характеристики турбулентной атмосферы. Описание метода

вычисления дисперсии

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

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

  1. Величина пульсации скорости в пределах объема, который занимает самолет обычных размеров, существенно не меняется.
  2. Пульсация скорости ветра является стационарным случайным процессом. Компоненты этой скоростиWx,Wy,Wz являются независимыми. Статистические характеристики пульсаций скорости ветра в поперечных направленияхWy,Wx одинаковы.
  3. Спектральные плотности компонентWx,Wy имеют следующие выражения:

  ;.

Где:

V-скорость самолета , м/с;  масштаб турбулентности ,м;

- частота порыва , 1/с;  - среднее квадратическое отклонение пульсации скорости ветра, м/с .

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

где,.- передаточные функции, с учетом замены .

Вычисление дисперсии сводится к вычислению несобственного интеграла

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

,

где , , А и В – полиномы с рациональными коэффициентами:

;

.

Необходимо, чтобы полином А(p) имел все нули в левой полуплоскости, а полиномB(p) имел все нули в левой полуплоскости и может быть на мнимой оси. Кроме того, степень полиномаB(p) должна быть по крайней мере, на единицу меньше, чем степень полинома А(p).

Если эти требования к полиномам выполнены, то дисперсия может быть вычислена по рекуррентному соотношению:

   ;k=1,2,…n  с начальным условием =0,

Здесь .

Параметры - коэффициенты полиномовAk(p) иBk(p),  степени которых не превосходятn;

  ;

  ;

  ;

  ;

;

  ;

.

Представление (факторизацию) числителя и знаменателя подынтегрального выражения в видеB(p)B(-p)иA(p)B(-p)можно получить следующим образом:

Квадрат модуля комплексного числа можно представить как

,

;

;

Где:

где:

;

Таблица коэффициентов:

;

;

;

;

;

;

;

.

Используемые табличные данные, характеризующие турбулентность:

H,км

/c

, м/c

L,км

1

1.65

1.36

0.624

2

1.65

1.43

0.831

4

2.04

1.68

0.972

6

2.13

1.69

1.01

8

2.15

1.69

0.98

10

2.23

1.73

1.1

12

2.47

1.79

1.54

Производную нормальной перегрузки по углу атаки найдем как:

, ,  , , ,

, , ,

,

Программная реализация.

Учитывая все вышеизложенное, разработаем программу вMatLab. Ее исходный код с комментариями приведен ниже:

Основная программа:

% Самолет "Ягуар", фронтовой истребитель

% 01-414,Лавренюк П.А.

clear all;

%ИСХОДНЫЕ ДАННЫЕ:

% Площадь несущей поверхности, м^2

S = 25;

% САХ, м

   b_a = 2.8;

% Момент инерции относительно осиZ, кг * м^2:

I_z = 2.05*10^5;

% Относительное положение центра тяжести:

X_t_otn = 0.3;

% Взлетная масса самолета, кг

m0 = 15000;

% Ускорение свободного падения, м / с^2

g = 9.81;

% Расстояние от центра тяжести до кабины, м

Lk = 5

% Атмосферное давление как функция высоты, Па

P_H = [79490 61660 47210 35650 26290 19390];

% Скорость звука как функция высоты, м / с

a_H = [332.7 324.7 316.6 308.2 299.2 295.2];

% Плотность воздуха как функция высоты, кг / м^3

ro_H = [1.01 0.819 0.66 0.526 0.41 0.31];

% Массив производной коэф. подъемной силы как функции М.

C_y_alpha = [3.6 3.9 4 4.1 4.2 4.3 4.5 4.7 4.8];

% Относительное положение фокуса как функции М:

X_f_otn = [0.36 0.4 0.43 0.425 0.42 0.45 0.5 0.56 0.6];

% Частная производная коэф. момента тангажа по отклонению руля высоты:

m_z_delta_v = [-1.126  -1.247 -1.249 -1.248 -1.247 -1.248 -1.251 -1.4 -1.35];

% Частная производная коэф. момента тангажа по угловой скорости:

m_z_omega_z = [-3.3 -3.3 -3.3 -3.3 -3.3 -3.4 -3.6 -4 -3.75];

% Частная производная коэф. момента тангажа по производной угла атаки:

m_z_alpha_t = [-1.06 -1.04 -1.02 -0.98 -0.92 -0.87 -0.87 -0.62 -0.55];

% СЛЕДУЮЩИЕ ПАРАМЕТРЫ В ПРОЦЕССЕ РАБОТЫ ПРОГРАММЫ МОГУТ БЫТЬ ПЕРЕГРУЖЕНЫ

% Масштаб турбулентности, м:

L = 1000;

% Среднее квадратическое отклонение пульсаций горизонтальной

% и вертикальной составляющих скоростей ветра:

Sigma_X = 5;

       Sigma_Y = 5;

% Коэффициент усиления по перегрузке:

K_n_y = 0.1:0.1:1;

% Коэффициент усиления по угловой скорости:

K_omega_z = 0.5:0.5:5;

% Числа М полета:

M = 0.5:0.1:1.3;

% Высоты полета, м:

H = 2000:2000:12000;

% Ввод пользователем аргумента(-тов) дисперсии:

disp('Введите порядковый номер зависимости, которую следует построить:');

disp('1 - D(K_n_y, K_omega_z), 2 - D(H), 3 - D(L), 4 - D(M), 5 - D(Sigma_Y),');

  disp('6 - Расчет контрольной точки для оценки гор. и верт. составляющих дисперсии,');

disp('7 - Расчет дисперсии приK_n_y = 0,K_omega_z = 0,');

  arg_of_disp = input('');

% Проверка допустимости введенного значения:

if or(arg_of_disp < 0, arg_of_disp > 7)

  disp('Недопустимый ввод! Программа будет завершена.');

  clear all;

  return;

end;

%РАСЧЕТ МАССИВОВ СКОРОСТИ И ПРОИЗВОДНЫХ МОМЕНТА ТАНГАЖА И НОРМ. ПЕРЕГРУЗКИ ПО УГЛУ АТАКИ

% Внешний цикл по высоте -i, внутренний цикл по М -j:

for i=1:1:6

  for j=1:1:9

     V(i, j) = M(j) * a_H(i);

     n_y_alpha(i, j) = (C_y_alpha(j) * ro_H(i) * V(i, j)^2 * S) / (2 * m0 * g);

     M_z_alpha(i, j) = ...

        (X_t_otn - X_f_otn(j)) * C_y_alpha(j) * ro_H(i) * V(i, j)^2 * S * b_a / (2 * I_z);

     M_z_delta_v(i, j) = ...

        m_z_delta_v(j) * ro_H(i) * V(i, j)^2 * S * b_a / (2 * I_z);

     M_z_omega_z(i, j) = ...

        (m_z_omega_z(j) * ro_H(i) * V(i, j)^2 * S * b_a^2) / (2 * I_z * V(i, j));

     M_z_alpha_t(i, j) = ...

        (m_z_alpha_t(j) * ro_H(i) * V(i, j)^2 * S * b_a^2) / (2 * I_z * V(i, j));

  end;

end;

% Вывод рассчитанных выше массивов в командное окноMatLab для проверки

V

n_y_alpha

M_z_alpha

M_z_delta_v

M_z_omega_z

M_z_alpha_t

% ВЫБОР АЛГОРИТМА РАСЧЕТА В ЗАВИСИМОСТИ ОТ ВВЕДЕННОГО ПОЛЬЗОВАТЕЛЕМ ЗНАЧЕНИЯ:

switch arg_of_disp

     % D(K_n_y, K_omega_z)

  case 1

i = 2; % Фиксация высоты - Н = 4000 м

j = 3; % Фиксация скорости - М = 0.7

forrow = 1:1:10 % Цикл по элементам массиваK_n_y

forcolumn = 1:1:10 % Цикл по элементам массиваK_omega_z

              % Вызов функции расчетов полиномовB_y(p) и А_y(p)

[B_y, A_y] = y_coeffs(i, j, row, column, Sigma_Y, L, V, n_y_alpha, ...

                 M_z_alpha, M_z_delta_v, M_z_omega_z, M_z_alpha_t, K_omega_z, K_n_y, g);

% Установка максимального порядка полиномов

              N = 4;

              % Вызов функции расчета дисперсии:

              result(row, column) = dispersion(A_y, B_y, N);

           end;

        end;

        % Построение поверхности D(K_n_y, K_omega_z):

        [K_n_y, K_omega_z] = meshgrid(0.1:0.1:1, 0.5:0.5:5);

           surf(K_n_y, K_omega_z, result);

           hold on;

           title('D(K\_n\_y, K\_omega\_z), L = 1000, M = 0.7, H = 4000, Sigma\_Y = 5');

           xlabel('K\_n\_y');

           ylabel('K\_omega\_z');

           zlabel('D(K\_n\_y, K\_omega\_z)');

           hold off;

% D(H)

case 2

j = 3; % Фиксация скорости - М = 0.7

row = 9; column = 9; % Фиксация K_n_y и K_omega_z

        fori=1:1:6% Цикл по Н

              % Вызов функции расчетов полиномовB_y(p) и А_y(p)

[B_y, A_y] = y_coeffs(i, j, row, column, Sigma_Y, L, V, n_y_alpha, ...

                 M_z_alpha, M_z_delta_v, M_z_omega_z, M_z_alpha_t, K_omega_z, K_n_y, g);

% Установка максимального порядка полиномов

              N = 4;

              % Вызов функции расчета дисперсии:

result(i) = dispersion(A_y, B_y, N);

         end;

         % Построение графикаD(H):

      plot(H, result, '-p');

           hold on;

           title('D(H), K\_n\_y = 1, K\_omega\_z = 5, L = 1000, M = 0.7, Sigma\_Y = 5');

           xlabel('H');

           ylabel('D(H)');

           grid on;

           hold off;

     % D(L)

     case 3

        i = 2; % Фиксация высоты - Н = 4000 м

j = 3; % Фиксация скорости - М = 0.7

row = 9; column = 9; % Фиксация K_n_y и K_omega_z

% ПерегрузкаL:

L = 0:1:1000;

L_array_index = 0; % Индекс элементов массиваL

for L_array_index = 1:1:1001 % Цикл по Н

% Вызов функции расчетов полиномовB_y(p) и А_y(p)

[B_y, A_y] = y_coeffs(i, j, row, column, Sigma_Y, L(L_array_index), ...

                    V, n_y_alpha, M_z_alpha, M_z_delta_v, M_z_omega_z, M_z_alpha_t, ...

                    K_omega_z, K_n_y, g);

% Установка максимального порядка полиномов

              N = 4;

              % Вызов функции расчета дисперсии:

result(L_array_index) = dispersion(A_y, B_y, N);

         end;

         % Построение графикаD(L):

         plot(L, result);

           hold on;

           title('D(L), K\_n\_y = 1, K\_omega\_z = 5, M = 0.7, H = 4000, Sigma\_Y = 5');

           xlabel('L');

           ylabel('D(L)');

           grid on;

           hold off;

     % D(M)

case 4

        i = 2; % Фиксация высоты - Н = 4000 м

row = 9; column = 9; % Фиксация K_n_y и K_omega_z

        forj=1:1:9 % Цикл поM

              % Вызов функции расчетов полиномовB_y(p) и А_y(p)

[B_y, A_y] = y_coeffs(i, j, row, column, Sigma_Y, L, V, n_y_alpha, ...

                 M_z_alpha, M_z_delta_v, M_z_omega_z, M_z_alpha_t, K_omega_z, K_n_y, g);

% Установка максимального порядка полиномов

              N = 4;

              % Вызов функции расчета дисперсии:

result(j) = dispersion(A_y, B_y, N);

         end;

         % Построение графикаD(M):

         plot(M, result, '-p');

           hold on;

           title('D(M), K\_n\_y = 1, K\_omega\_z = 5, L = 1000, H = 4000, Sigma\_Y = 5');

           xlabel('M');

           ylabel('D(M)');

           grid on;

           hold off;

     % D(Sigma_Y)

case 5

i = 2; % Фиксация высоты - Н = 4000 м

j = 3; % Фиксация скорости - М = 0.7

row = 9; column = 9; % Фиксация K_n_y и K_omega_z

% ПерегрузкаSigma_Y:

Sigma_Y = 0.5:0.5:5;

Sigma_Y_array_index = 0; % Индекс элементов массиваSigma_Y

for Sigma_Y_array_index = 1:1:10 % Цикл по Sigma_Y

% Вызов функции расчетов полиномовB_y(p) и А_y(p)

[B_y, A_y] = y_coeffs(i, j, row, column, Sigma_Y(Sigma_Y_array_index), ...

                    L, V, n_y_alpha, M_z_alpha, M_z_delta_v, M_z_omega_z, ...

                    M_z_alpha_t, K_omega_z, K_n_y, g);

% Установка максимального порядка полиномов

              N = 4;

              % Вызов функции расчета дисперсии:

result(Sigma_Y_array_index) = dispersion(A_y, B_y, N);

         end;

         % Построение графикаD(Sigma_Y):

         plot(Sigma_Y, result, '-p');

           hold on;

           title('D(Sigma\_Y), K\_n\_y = 1, K\_omega\_z = 5, L = 1000, H = 4000, M = 0.7');

           xlabel('Sigma\_Y');

           ylabel('D(Sigma\_Y)');

           grid on;

           hold off;

% Расчет контрольной точки

case 6

i = 2; % Фиксация высоты - Н = 4000 м

j = 3; % Фиксация скорости - М = 0.7

row = 9;column = 9; % Индексы элементов массивовK_n_y иK_omega_z

        % Вызов функции расчетов полиномовB_y(p) и А_y(p)

     [B_y, A_y] = y_coeffs(i, j, row, column, Sigma_Y, L, V, n_y_alpha, ...

           M_z_alpha, M_z_delta_v, M_z_omega_z, M_z_alpha_t, K_omega_z, K_n_y, g);

% Установка максимального порядка полиномов дляB_y(p) и А_y(p)

        N = 4;

        % Вызов функции расчета дисперсии:

disp('Дисперсия вертикальной составляющей: ');

        result_y = dispersion(A_y, B_y, N)

% Вызов функции расчетов полиномовB_x(p) и А_x(p)

     [B_x, A_x] = x_coeffs(i, j, row, column, Sigma_X, L, V, n_y_alpha, ...

           M_z_alpha, M_z_delta_v, M_z_omega_z, M_z_alpha_t, K_omega_z, K_n_y, g);

% Установка максимального порядка полиномов дляB_x(p) и А_x(p)

        N = 3;

        % Вызов функции расчета дисперсии:

disp('Дисперсия горизонтальной составляющей: ');

        result_x = dispersion(A_x, B_x, N)

% Расчет дисперсии приK_n_y = 0,K_omega_z = 0

case 7

i = 2; % Фиксация высоты - Н = 4000 м

j = 3; % Фиксация скорости - М = 0.7

row = 1;column = 1; % Индексы элементов массивовK_n_y иK_omega_z

        % ПерегрузкаK_n_y иK_omega_z:

K_n_y = 0;

        K_omega_z = 0;

% Вызов функции расчетов полиномовB_y(p) и А_y(p)

     [B_y, A_y] = y_coeffs(i, j, row, column, Sigma_Y, L, V, n_y_alpha, ...

           M_z_alpha, M_z_delta_v, M_z_omega_z, M_z_alpha_t, K_omega_z, K_n_y, g);

% Установка максимального порядка полиномов

        N = 4;

        % Вызов функции расчета дисперсии и вывод результатов расчетов вMatLab:

disp('Дисперсия приK_n_y = 0,K_omega_z = 0 равна: ');

result = dispersion(A_y, B_y, N)

     % Перегрузка K_n_y и K_omega_z:

          K_n_y = 0.1;

         K_omega_z = 0.5;

        % Вызов функции расчетов полиномовB_y(p) и А_y(p)

     [B_y, A_y] = y_coeffs(i, j, row, column, Sigma_Y, L, V, n_y_alpha, ...

           M_z_alpha, M_z_delta_v, M_z_omega_z, M_z_alpha_t, K_omega_z, K_n_y, g);

% Установка максимального порядка полиномов

        N = 4;

        % Вызов функции расчета дисперсии и вывод результатов расчетов вMatLab:

        disp('Дисперсия при малыхK_n_y = 0.1,K_omega_z = 0.5 равна: ');

result = dispersion(A_y, B_y, N)

        

Функция расчета коэффициентов полиномовB_x(p) и А_x(p)

function [B_x, A_x] = x_coeffs(i, j, row, column, Sigma_X, L, V, n_y_alpha, ...

  M_z_alpha, M_z_delta_v, M_z_omega_z, M_z_alpha_t, K_omega_z, K_n_y, g)

% Расчет коэффициентов полиномаB_x(p):

  b_x_0(row, column) = Sigma_X * sqrt(L / (pi * V(i, j))) * ...

     M_z_alpha_t(i, j) / n_y_alpha(i, j);

     b_x_1(row, column) = Sigma_X * sqrt(L / (pi * V(i, j))) * ...

        (M_z_alpha(i, j) / n_y_alpha(i, j) + ...

        M_z_delta_v(i, j) * K_n_y(row));

% Расчет коэффициентов полиномаA_x(p):

  a_x_0(row, column) = L / (V(i, j) * n_y_alpha(i, j));

     a_x_1(row, column) = 1 / n_y_alpha(i, j) + L * (g * n_y_alpha(i, j) - ...

        V(i, j) * M_z_omega_z(i, j) - V(i, j) * M_z_delta_v(i, j) * ...

        K_omega_z(column) - V(i, j) * M_z_alpha_t(i, j)) / ...

        (V(i, j)^2 * n_y_alpha(i, j));

     a_x_2(row, column) = (g * n_y_alpha(i, j) - V(i, j) * ...

        M_z_omega_z(i, j) - V(i, j) * M_z_delta_v(i, j) * K_omega_z(column) ...

        - V(i, j) * M_z_alpha_t(i, j)) / (V(i, j) * n_y_alpha(i, j)) + L * ...

        (-g * n_y_alpha(i, j) * M_z_omega_z(i, j) - g * n_y_alpha(i, j) * ...

        M_z_delta_v(i, j) * K_omega_z(column) - M_z_alpha(i, j) * V(i, j) - ...

        M_z_delta_v(i, j) * K_n_y(row) * n_y_alpha(i, j) * V(i, j)) / ...

        (V(i, j) * n_y_alpha(i, j));

     a_x_3(row, column) = (-g * n_y_alpha(i, j) * M_z_omega_z(i, j) - ...

        g * n_y_alpha(i, j) * M_z_delta_v(i, j) * K_omega_z(column) - ...

     M_z_alpha(i, j) * V(i, j) - M_z_delta_v(i, j) * K_n_y(row) * ...

        n_y_alpha(i, j) * V(i, j)) / (V(i, j) * n_y_alpha(i, j));

B_x = [b_x_0(row, column) b_x_1(row, column) 0];

A_x = [a_x_0(row, column) a_x_1(row, column) a_x_2(row, column) ...

     a_x_3(row, column)];

Функция расчета коэффициентов полиномов B_y(p) и А_y(p)

function [B_y, A_y] = y_coeffs(i, j, row, column, Sigma_Y, L, V, n_y_alpha, ...

M_z_alpha, M_z_delta_v, M_z_omega_z, M_z_alpha_t, K_omega_z, K_n_y, g)

% Расчет коэффициентов полинома B_y(p):

b_y_0(row, column) = (Sigma_Y * sqrt((6 * L) / (pi * V(i, j))) * L) ...

     / (2 * V(i, j)^2);

     b_y_1(row, column) = (Sigma_Y /(2 * V(i, j))) * sqrt((2 * L) / ...

        (pi * V(i, j))) * (1 - sqrt(3) * L * (M_z_omega_z(i, j) + ...

        M_z_delta_v(i, j) * K_omega_z(column)) / V(i, j));

b_y_2(row, column) = (Sigma_Y/2) * sqrt((2 * L) / (pi * V(i, j))) * ...

  (M_z_omega_z(i, j) + M_z_delta_v(i, j) * K_omega_z(column)) / V(i, j);

% Расчет коэффициентов полинома A_y(p):

  a_y_0(row, column) = L^2 / (V(i, j)^2 * n_y_alpha(i, j));

     a_y_1(row, column) = (2 * L) / (V(i, j) * n_y_alpha(i, j)) + (L^2 * ...

        (g * n_y_alpha(i, j) - V(i, j) * M_z_omega_z(i, j) - V(i, j) * ...

        M_z_delta_v(i, j) * K_omega_z(column) - ...

        M_z_alpha_t(i, j) * V(i, j))) / (V(i, j)^3 * n_y_alpha(i, j));

     a_y_2(row, column) = 1 / n_y_alpha(i, j) + 2 * L * (g * ...

        n_y_alpha(i, j) - V(i, j) * M_z_omega_z(i, j) - V(i, j) * ...

        M_z_delta_v(i, j) * K_omega_z(column) - ...

        M_z_alpha_t(i, j) * V(i, j)) / (V(i, j)^2 * n_y_alpha(i, j));

     a_y_3(row, column) = (g * n_y_alpha(i, j) - ...

        V(i, j) * M_z_omega_z(i, j) - V(i, j) * M_z_delta_v(i, j) * ...

        K_omega_z(column) - M_z_alpha_t(i, j) * V(i, j)) / (V(i, j) * ...

        n_y_alpha(i, j)) + 2 * L * (-g * n_y_alpha(i, j) * ...

        M_z_omega_z(i, j) - g * n_y_alpha(i, j) * M_z_delta_v(i, j) * ...

        K_omega_z(column) - M_z_alpha(i, j) * V(i, j) - M_z_delta_v(i, j) ...

        * K_n_y(row) * n_y_alpha(i, j) *  V(i, j)) / ...

        (V(i, j)^2 * n_y_alpha(i, j));

     a_y_4(row, column) = (-g * n_y_alpha(i, j) * M_z_omega_z(i, j) - g * ...

        n_y_alpha(i, j) * M_z_delta_v(i, j) * K_omega_z(column) - ...

        M_z_alpha(i, j) * V(i, j) - M_z_delta_v(i, j) * ...

        K_n_y(row) * n_y_alpha(i, j) *  V(i, j)) / (V(i, j) * ...

        n_y_alpha(i, j));

B_y = [b_y_0(row, column) b_y_1(row, column) b_y_2(row, column) 0];

A_y = [a_y_0(row, column) a_y_1(row, column) a_y_2(row, column) ...

     a_y_3(row, column) a_y_4(row, column)];

Функция расчета дисперсии:

function s=v(an,bn,n)

a=an;b=bn;

pi=3.14;

c=0;

ier=0;

if a(1)<=0

   s = NaN;

ier=1

   return

end

for k=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

       s=1000000

       ier=1

       return

   end

end

s=pi*c ;

Результаты расчетов.

  1. Зависимость дисперсии от коэффициентов и . В качестве опорного режима полета здесь и далее выбираем: М = 0.7, Н = 4000 м.

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

Из графика видно, что с ростом  дисперсия перегрузки увеличивается, а

при возрастании  -  убывает.

  1. Зависимость дисперсии от высоты полета.

Мы видим, что с ростом высоты полета дисперсия перегрузки убывает.

  1. Зависимость дисперсии от числа М полета:

С увеличением числа М дисперсия перегрузки монотонно возрастает.

  1. Зависимость дисперсии от масштаба турбулентности:

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

  1. Зависимость дисперсии от среднеквадратического отклонения вертикальной составляющей скорости ветра:

Дисперсия  перегрузки возрастает с увеличением.

  1. Расчет контрольной точки.

Данный  расчет проводится для того, чтобы оценить вклад горизонтальной и вертикальной составляющих дисперсии в формуле (4.3). Программа производит расчет указанных составляющих при следующих параметрах полета: М = 0.7, Н = 4000 м,L = 1000 м, ,,.

В результате расчета получим:

Дисперсия вертикальной составляющей:

result_y = 0.1101

Дисперсия горизонтальной составляющей:

result_x = 0.0606

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

  1. Расчет дисперсии на опорном режиме при, .

Здесь произведем расчет дважды: при,  и при , .

В результате расчета получим:

Дисперсия при K_n_y = 0, K_omega_z = 0 равна:

result =

  NaN

Дисперсия при малых K_n_y = 0.1, K_omega_z = 0.5 равна:

result =

   0.4491

Вывод: с уменьшением  и дисперсия стремится к нулю.

Список используемой литературы.

  1. Аэромеханика самолета. Под редакцией Бочкарева А. Ф.  и Андреевского В. В.
  2. Остославский И. В. Аэродинамика самолета.
  3. Гуськов Ю. П., Загайнов Г. И. Управление полетом самолетов.
  4. Овчаренко В. Н, Павлов К. А. Методические указания к курсовой работе по теме «Статистическая динамика».
  5. Маркин Н. Н. Курс лекций по теории автоматического управления.
  6. Оглоблин А. В. Курс лекций по динамике полета.
  7. Овчаренко В. Н. Курс лекций по теории вероятностей и статистической динамике.


 

А также другие работы, которые могут Вас заинтересовать

41096. НЕОБХОДИМОСТЬ ДЕНЕГ, ИХ ВОЗНИКНОВЕНИЕ И СУЩНОСТЬ 656.5 KB
  Деньги возникают при определенных условиях осуществления производства и экономических отношений в обществе и способствуют дальнейшему их развитию.
41097. СИСТЕМА БЕЗНАЛИЧНЫХ РАСЧЕТОВ 627.73 KB
  Сущность принципы организации и значение безналичных расчетов. Аккредитивная форма расчетов ее сущность и сфера применения. Денежные средства на расчетных и других аналогичных счетах в банках отражаются посредством записи остатков оборотов по лицевым счетам вследствие безналичных расчетов.
41098. Коммерческие банки. Сущность и организационная основа деятельности коммерческих банков 103.42 KB
  Принципы деятельности коммерческого банка. Функции коммерческого банка Банки одно из центральных звеньев системы рыночных структур. Основное назначение банка –посредничество в перемещении денежных средств от кредиторов к заемщикам и от продавцов к покупателям. Наряду с банками перемещение денежных средств на рынках осуществляют и другие финансовые и кредитнофинансовые учреждения: инвестиционные фонды страховые компании брокерские дилерские фирмы и т.
41099. Управление заемным капиталом 1.25 MB
  Обеспечение своевременных расчетов по полученным кредитам На второй стадии анализа определяются основные формы привлечения заемных средств анализируются в динамике удельный вес сформированных финансового кредита товарного кредита и текущих обязательств по расчетам в общей сумме заемных средств используемых предприятием. Эти формы дифференцируются в разрезе финансового кредита; товарного коммерческого кредита; прочих форм. К числу важнейших из этих условий относятся; а срок предоставления кредита; б ставка процента за кредит;...
41101. ОСНОВЫ МЕЖДУНАРОДНЫХ ВАЛЮТНЫХ, РАСЧЕТНЫХ И КРЕДИТНЫХ ОТНОШЕНИЙ 370.33 KB
  С середины 70х годов используется паритет на базе валютной корзины; режим валютного курса фиксированный плавающий в определенных пределах; наличие или отсутствие валютных ограничений; регулирование международной валютной ликвидности страны; официальные золотые и валютные резервы стран; счета СДР ЭКЮ резервная позиция в МВФ; регламентация использования международных кредитных средств обращения и форм международных расчетов; режимы валютного рынка и рынка золота;статус национальных органов регулирующих валютные...
41102. Центральные и коммерческие банки, их функции и формы организации 292.94 KB
  Введение Термин банк происходит от итальянского слова банко что означает лавка скамья или конторка за которой менялы оказывали свои услуги. Во многих источниках дошедших до нас можно встретить данные о вавилонских банкирах принимавших процентные вклады и выдававших ссуды под письменные обязательства и под залог различных ценностей. Вавилонский банк принимал вклады платил по ним проценты выдавал ссуды и даже выпускал банковские билеты.
41103. Компьютерные мониторы на основе электронно-лучевой трубки 839 KB
  Сквозь металлическую маску или решетку они попадают на внутреннюю поверхность стеклянного экрана монитора которая покрыта разноцветными люминофорными точками. Поток электронов луч может отклоняться в вертикальной и горизонтальной плоскости что обеспечивает последовательное попадание его на все поле экрана. Чтобы электроны беспрепятственно достигали экрана из трубки откачивается воздух а между пушками и экраном создаётся высокое электрическое напряжение ускоряющее электроны.Это сделано для того чтобы электронный луч в центре экрана и...
41104. Цифровая бумага 1.87 MB
  Многоцветная полихромная электронная бумага Электронная бумага EDP В отличие от традиционных жидкокристаллических плоских дисплеев в которых используется просвет матрицы для формирования изображения электронная бумага формирует изображение в отраженном свете как обычная бумага и может показывать текст и графику неопределенно долго не потребляя при этом электричество и позволяя изменять изображение в дальнейшем.