Динамика полета воздушного судна

самолет вертолет аэродинамический

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

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

Цели (для самолёта):

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

Задачи (для самолёта):

1. Пересчитать аэродинамических характеристик профиля крыла самолёта, на крыло конечного размаха.

2. Определить минимальную тягу.

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

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

5. Решение системы алгебраических уравнений описывающих движение самолёта.

Цели (для вертолёта):

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

1. Построить 3D модель фюзеляжа вертолёта и лопастей несущего винта.

2. Построить сетку из конечных элементов для фюзеляжа и лопасти.

3. «Продуть» лопасть в программном пакете «ANSYS CFX».

4. Перенести поле скоростей от лопасти в модель фюзеляжа.

5. «Продуть» фюзеляж в программном пакете «ANSYS CFX».

6. С помощью выражений нормальной и продольной сил вычислить коэффициенты нормальной и продольной сил.

7. Вычислить коэффициенты подъёмной силы и силы лобового сопротивления из коэффициентов нормальной и продольной сил.

1. Самолёт

1.1 История создания самолёта, его массо-геометричекие и летно-технические характеристики

Прототип самолета Vickers Viastra Mk. I совершил свой первый полет 1 октября 1930 года. Это был десятиместный пассажирский самолет оснащенный тремя двигателями Armstrong Siddeley Lynx Major мощностью 270 л.с. После испытаний самолет значительно модернизировали, заменив три двигателя на два более мощных Bristol Jupiter XIF (520 л.с.) и расширив пассажирскую кабину еще на два места.

Viastra Mk. I это цельнометаллический моноплан, главной особенность которого является использование трёх двигателей Armstrong Siddeley Lynx Major мощностью 270 л.с. Два двигателя довешаны на крыльях, третий встроен в носовую часть самолёта.

Крылья самолёта прямоугольные и обшиты гофрированным металлом в них же находятся топливные баки, рассчитанные на полёт до 480 км. при крейсерской скорости. Основную нагрузку в крыле воспринимают лонжероны, которые на фланцах имеют дюралюминиевые полосы для компенсации местных напряжений. Центральные части крыльев соединены с фюзеляжем двумя трубами

В основе фюзеляжа лежит скелетный каркас, который подкреплён дюралюминиевой приклёпанной обшивкой для большей устойчивости. Каркас состоит из Т-образных лонжеронов. Пол кабины укрепляется поперечными опорами.

Пустым самолёт весит 3170 килограмм, а с полезной нагрузкой 4840 килограмм, он развивал крейсерскую скорость в 193 км/ч, что в те времена было большим достижением для самолёта с тремя двигателями.

Cамолет поставлялся австралийской авиакомпании Western Australian Airways в 1931 году. Последний из австралийских Viastra был списан в 1936 года. Вслед за первыми серийными самолетами, последовали несколько модификаций выполненных на заказ. Такими были Viastra Mk. VI с одним двигателем Jupiter XIF, и Viastra Mk. X, построенный специально по заказу Принца Уэльского, летавшего на нем до 1934 года.

Таблица 1 — Массо-геометричекие и летно-технические

Рисунок 1 — Геометрические характеристики самолёта

1.2 Аэродинамические характеристики профиля RAF-34

Профиль является не симметричным. Его модификации использовались в крыльях самолётов Ту-4, ТИ-28, в британском транспортном планере General Aircraft G.A.L.49 Hamilcar.

Рисунок 2 — Профиль Raf-34

Рисунок 4 — Коэффициенты подъёмной силы и лобового сопротивления профиля Рисунок 5 — Качество профиля Таблица 2 — Аэродинамические характеристики[2]

1.3 Аэродинамические характеристики крыла конечного размаха

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

При наличии подъёмной силы под крылом образуется зона повышенного давления, а над крылом — пониженного давления. Из-за разности давлений происходит перетекание воздушного потока через торцы крыла с нижней поверхности на верхнюю, что приводит к появлению дополнительного потока, направленного вдоль размаха крыла. Под крылом этот поток направлен к концам крыла, а над ним — к его средней части [3,4].

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

По теореме Н. Е. Жуковского индуктивное сопротивление определяется следующим выражение [3]:

где

— удлинение крыла;

b=3,87 м — хорда крыла;

l=21,33 м — размах крыла;

— коэффициент, учитывающий форму крыла.

Сравним коэффициенты лобового сопротивления с учётом индуктивного сопротивления и без него по рисунку.

Сравнение коэффициентов лобового сопротивления

1.4 Аэродинамика крыла конечного размаха с учётом механизации

Механизация крыла необходима для улучшения лётно-технических характеристик самолёта на режимах разбега, взлёта посадки и набора высоты. Механизация даёт возможности:

— уменьшить посадочную скорость;

— уменьшить длину разбега;

— улучшить устойчивость и управляемость на больших углах атаки.

На данном самолёте из механизации имеются щитки. Это управляемые поверхности, прилегающие к задней поверхности крыла. При их отклонении повышается коэффициент подъёмной силы в среднем на 50%, а коэффициент лобового сопротивления на 18%.

В итоге аэродинамические коэффициенты крыла имеют следующий вид:

Коэффициенты подъёмной силы и лобового сопротивления крыла с механизацией и без неё

Качество крыла с механизацией

1.5 Определение оптимальных параметров движения самолёта

Для определения оптимальных параметров движения самолёта воспользуемся методом тяг Жуковского[3,4].

Качество будет иметь максимально значение при наивыгоднейшем угле касательной к её графику. Тогда качество определится как

а наивыгоднейший угол атаки

.

Минимальную тягу и наивыгоднейшую скорость вычислим по формулам в несколько итераций:

Первая итерация:

Н,

м/с, Вторая итерация:

Н,

м/с, Третья итерация:

Н,

м/с, Номинальная тяга:

Н Максимальная тяга:

Н, Результаты приведены в таблице 3.

Таблица 3 — Расчётные параметры горизонтального полёта

1.6 Динамика самолёта

Определить движение самолёта — значит определить все его параметры полёта: положение самолёта, его скорость, высоту, ускорение, угловую скорость и т. п. в каждый момент времени в заданной системе координат.

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

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

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

1.7 Определение траекторных параметром движения самолёта

При определении траекторных параметров движения самолёта примем следующими допущениями:

1) угловых скоростей самолёта относительно всех его осей равны нулю,

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

3) механизация крыла всегда включена,

4) тяга всегда имеет максимальное значение.

Система уравнений движения центра масс самолёта по траектории

С учётом принятых допущений система дифференциальных уравнений движения центра масс самолёта по траектории в стартовой системе координат будет иметь вид:

где

— лобовое сопротивление самолёта в скоростной системе координат;

— подъёмная сила самолёта в скоростной системе координат;

м² — площадь крыла;

— скоростной напор;

— суммарная скорость;

— угол между горизонтом и вектором скорости V;

— угол атаки;

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

— функция от высоты;

м — радиус Земли;

м/с² — ускорение свободного падения;

— функция плотности от высоты;

кг/м³ — плотность воздуха на уровне моря (y=0 м);

— тяга;

— максимальная тяга;

Н — минимальная тяга;

— сила реакции опоры;

— сила трения;

— коэффициент трения для сухого бетона;

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

Метод конечных разностей

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

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

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

Из формулы (10), после отбрасывания третьего и последующих слагаемых, производная преобразуется в выражение:

которое является дискретным аналогом «вверх по потоку» первого порядка точности.

Система алгебраических уравнений

Используя выражение (11), преобразуем системы дифференциальных уравнений (10) в систему алгебраических уравнений:

где

— сила тяжести на i-м шаге решения,

— функция от высоты на i-м шаге решения,

— функция плотности от высоты на i-м шаге решения,

— тяга двигателей на i-м шаге решения,

— угол атаки самолёта на i-м шаге решения,

— сила трения на i-м шаге решения,

— силы реакции опоры на i-м шаге решения,

— угол между горизонтом и вектором скорости V на i-м шаге решения,

— суммарная скорость на i-м шаге решения,

— скоростной напор на i-м шаге решения,

— сила лобового сопротивления на i-м шаге решения,

— подъёмная сила на i-м шаге решения.

Для решения системы алгебраических уравнений (14) был написан скрипт в «Matlab» (см. приложение 1). Решение производилось с шагом по времени dt=0.1c, количество итераций составило 80 000. Время расчёта — 20 секунд.

Время разбега определим по следующему условию:

.

Таблица 4 — Результаты расчёта

Графики зависимостей аэродинамических коэффициентов, скорости, угла атаки, силы трения ускорения, плотности от времени и траектория представлены на рисунках 8−16.

Траектория движения самолёта График зависимости скорости от времени График зависимости угла атаки от времени График зависимости коэффициента лобового сопротивления от времени График зависимости коэффициента подъёмной силы от времени График зависимости силы трения от времени при разбеге самолёта График зависимости ускорения от времени в течение 100 с полёта График изменения плотности по времени График изменения ускорения свободного падения по времени

1.8 Определение параметров движения самолёта с учётом вращения вкруг центра масс

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

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

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

где

— статический момент от стабилизатора,

— динамический момент от стабилизатора, возникающий из-за приращения к углу атаки (см. рисунок 18) и оказывающее демпфирующее воздействие при вращении самолёта,

— статический момент от руля высоты,

— коэффициент нормальной силы статического момента на рулях высоты,

— угол атаки руля высоты (рисунок 17),

— угол отклонения руля высоты от начально положения,

— коэффициент нормальной силы статического момента на стабилизаторе,

— коэффициент нормальной силы динамического момента на стабилизаторе,

— приращение к углу атаки на стабилизаторе (рисунок 18),

м — расстояние от центра масс самолёта до центра давления хвостового оперения,

м² — площадь руля высоты,

м² — площадь стабилизатора,

кг/м³ — момент инерции самолёта относительно оси OZ,

— угловая скорость вращения самолёта вокруг его центра масс в плоскости XY.

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

Метод Рунге-Кутта

Метод Рунге — Кутта непосредственно рассчитан на интегрирование дифференциальных уравнений вида:

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

Тогда значение функции в точке i+1 будет вычисляться как:

Дискретный аналог

Для решения дифференциального уравнения (16) «заморозим» величину t поскольку она не присутствует в выражении в явном виде, тогда функция (19) примет вид:

В итоге аналогично (23) получаем угловую скорость вращения самолёта вокруг центра масс в плоскости YX:

а угол тангажа:

Функцию, аналогичную (19), получим подстановкой (16) в (17):

где

— скоростной напор,

— аэродинамический коэффициент нормальной силы на рулях высоты на i-м шаге решения,

— угол атаки на рулях высоты на i-м шаге решения, Подставляя выражение (27) в (20) получим:

Подставляя выражение (27) в (21), с учётом приращений, получим:

где

— аэродинамический коэффициент нормальной силы на стабилизаторе на i-м шаге решения с приращением,

— угол атаки на стабилизаторе на i-м шаге решения с приращением.

Подставляя выражение (27) в (22), с учётом приращений, получим:

где

— аэродинамический коэффициент нормальной силы на стабилизаторе на i-м шаге решения с приращением,

— угол атаки на стабилизаторе на i-м шаге решения с приращением.

Подставляя (28) — (30) в (25), получим выражение угловой скорости для численного решения.

Тогда система дифференциальных уравнений (15) c учётом полученной ранее в п. 4.1.3. системой (14) преобразуется в систему алгебраических уравнений:

Параметры вращательного движения самолёта полученные при расчёте на ЭВМ

Для решения системы алгебраических уравнений (31) был написан скрипт в «Matlab» (см. приложение 2). Решение производилось с шагом по времени dt=0.01c, количество итераций составило 10 000. Время расчёта — 5 секунд. Результаты сведём в таблицу 4.

Таблица 5 — Результаты расчёта

Графики зависимостей угла тангажа, статического и динамического моментов стабилизатора, статического момента руля высоты от времени представлены на рисунках 19−24.

График зависимости угла атаки на руле высоты от времени в интервале от 4 до 16 с График зависимости динамического момента от стабилизатора от времени в интервале от 4 до 16 с График зависимости статического момента от стабилизатора от времени в интервале от 4 до 16 с График зависимости статического момента от руля высоты от времени в интервале от 4 до 16 с График зависимости угла тангажа от времени самолет вертолет аэродинамический График зависимости угла тангажа от времени в интервале от 20 до 21 с

2. Вертолёт

2.1 Балансировка вертолёта

Продольная балансировка вертолёта

Условием продольной балансировки вертолёта является равенство нулю суммы действующих на него продольных моментов при соблюдении равновесия проекций сил на оси OX и OY.

Силы, действующие на вертолёт в продольной плоскости Сумма продольных моментов складывается из моментов несущего винта, фюзеляжа, стабилизатора (рисунок 25):

Рассмотрим каждое слагаемое в отдельности.

Продольный момент несущего винта:

где

Т — тяга несущего винта,

Н — продольная сила несущего винта,

х_т — продольная центровка вертолёта,

у_т — вертикальная центровка вертолёта.

Момент аэродинамических сил фюзеляжа при поступательном полёте:

где

m_zф — коэффициент момента фюзеляжа,

S_усл — условная площадь, к которой отнесён момент фюзеляжа при определении m_zф,

L_усл — условный размер, к которому отнесён момент фюзеляжа при определении m_zф.

Продольный момент стабилизатора:

где

с_уст — коэффициент подъёмной силы стабилизатора,

S_ст — площадь стабилизатора,

L_ст — расстояние между центром тяжести вертолёта и центром давления стабилизатора.

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

Рассмотрим каждое слагаемое в отдельности.

G — силы тяжести,

v — угол тангажа,

ц_ст — установочный угол стабилизатора,

б — угол атаки.

Подъёмная сила фюзеляжа:

где

i — номер участка фюзеляжа,

n — число участков,

— площадь i — го участка в плане.

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

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

2.2 Боковая балансировка

Условиями боковой балансировки вертолёта являются равновесие моментов относительно вертикальной и продольной осей и равновесие проекций сил на поперечную ось:

Силы и моменты, действующие на вертолёт в п лоскости ZOY

Сумма моментов состоит из следующих слагаемых:

Силы и моменты, действующие на вертолёт в плоскости XOZ

Сумма моментов :

где

 — реактивные моменты несущего винта,

— момент, создаваемые рулевым винтом, вокруг ось OY,

— момент, создаваемые рулевым винтом, вокруг ось OX,

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

а сумма всех сил, действующих на вертолёт на ось OZ:

где

— боковая сила несущего винта,

— тяга рулевого винта,

— боковая сила планера,

— площадь боковой проекции фюзеляжа.

В итоге система уравнений боковой балансировки в связанной система координат будет иметь вид:

2.3 Аэродинамические характеристики вертолёта

Расчёт вертолёта в «ANSYS CFX»

3D-модель вертолёта, импортированную в «ANSYS ICEM», разбиваем за плоскости и разбиваем на конечные элементы так, чтобы их количество не превышало 1 млн.

Пример геометрии фюзеляжа вертолёта представлен на рисунке 27.

Геометрия фюзеляжа Пример сетки конечных элементов фюзеляжа вертолёта представлен на рисунке 28.

Разбиение фюзеляжа на конечные элементы Далее импортируем сетку в «ANSYS CFX», и назначаем граничные условия:

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

— стандартные параметры атмосферы на остальных плоскостях.

Граничные условия представлены на рисунке 29.

Граничные условия В раздел «expressions» вводим выражения для вычисления C_n и C_R:

.

Распределение давления по фюзеляжу Линии тока в плоскости YOZ

Линии тока в плоскости XOZ

Затем, изменяя угол атаки б перезапускаем расчёт. Результаты приведём в таблице 6.

Таблица 6 — Аэродинамические характеристики фюзеляжа

График зависимости аэродинамических коэффициентов от угла атаки

1 Журнал «Flight» — 26 сентября, 1930.

2 Кравец А. С. Характеристики авиационных профилей — М. «Государственное издательство оборонной промышленности», 1939 — 213 с.

3 Бочкарёв А. Ф. Аэромеханика самолёта/ Андреевский В. В. — М. «Машиностроение», 1985 — 361 с.

4 Байдаков В. Б. Аэродинамика и динамика полёта летательных аппаратов/ Клумов А. С. — М. «Машиностроение», 1979 — 342 с.

5 Флетчер К. Вычислительные методы в динамике жидкостей ч. 1 — М. «Мир», 1991 — 505 с.

6 Н. С. Бахвалов Численные методы/ Н. П. Жидков, Г. М. Кобельников — М. 1987 — 636 с.

7 С. В. Богословский Динамика полёта летательных аппаратов: Учеб. пособие/ А. Д. Дорофеев — СПб, СПбГУАП 2002 — 64 с.

8 А. С. Браверман Балансировки одновинтового вертолёта/ Д. М. Перлштейн, С. В. Лаписова — М. «Машиностроение», 1975 — 176 с.

9 А. С. Браверман Динамика вертолёта. Предельные режимы полёта/ А. П. Вайнтруб — М. «Машиностроение», 1988 — 275 с.

Приложение 1

%Начальные условия

m=4840;%масса самолёта

p0=1.25;%плотность воздуха при y=0

f=0.02;%коэффициент трения сухого бетона

dt=0.1;%шаг интегрирования

g0=9.8;%ускорение свободного падения при y=0

S=69.2;%площадь крыла самолёта

Pkr=3458.92;%минимальная тяга

PF=Pkr/(0.7*0.8);%максимальная тяга

n=80 000;%количество итераций

i=1;%счётчик

t (i)=0;%начальное время

y (i)=0;%начальная высота

x (i)=0;%начальное положение самолёта

teta (i)=0;%угол между вектором скорости и горизонтом

Phi (i)=8*pi/180;%угол тангажа

V0=0;%начальная скорость

%Решение

%Ускорение свободного падения

g (i)=g0*6 300 000² /(6 300 000+y (i))^2;

%Проекции начальной скорости на оси OX и OY

Vx (i)=V0*cos (teta (i));

Vy (i)=V0*sin (teta (i));

V (i)=((Vx (i)^2)+Vy (i)^2)^0.5;

%Углы

Alpha (i)=Phi (i) — teta (i);

%Плотность

po (i)=p0*(1-y (i)/44 300)^5.236;

%Тяга

P (i)=PF*po (i)/p0;

%Апроксимация аэродинамических коэффициентов

cyPoints=[-0.4 -0.84 -1.02 -0.63 0.06 0.63 1.02 0.84 0.4];

cxPoints=[0.015 0.04 0.11 0.2 0.32];

AlphaPoints=[-24 -20 -16 -8 0 8 16 20 24];

AlphaPoints2=[0 8 16 20 24];

p1=polyfit (AlphaPoints, cyPoints, 3);

p2=polyfit (AlphaPoints2, cxPoints, 3);

%Аэродинамические коэффициенты

Cy (i)=polyval (p1, Alpha (i)*180/pi)+0.6*polyval (p1, Alpha (i)*180/pi);

Cxi (i)=0.056*Cy (i)^2;

Cx (i)=polyval (p2, Alpha (i)*180/pi)+ Cxi (i);

%Аэродинамические силы

X (i)=po (i)*(V (i)^2)*Cx (i)*S/2;

Y (i)=po (i)*(V (i)^2)*Cy (i)*S/2;

%Сила реакции опоры

N (i)=m*g (i) — Y (i)*cos (teta (i))+X (i)*sin (teta (i)) — P (i)*sin (Phi (i));

%Сила трения

Ftr (i)=(N (i) — Y (i)*cos (teta (i)))*f;

PhiGrad (i)=Phi (i)*180/pi;

AlphaGrad (i)=Alpha (i)*180/pi;

b=1;

ay (i)=(P (i)*sin (Phi (i)) — m*g (i) — X (i)*sin (teta (i))+Y (i)*cos (teta (i)))/m;

ax (i)=(P (i)*cos (Phi (i)) — X (i)*cos (teta (i)) — Y (i)*sin (teta (i)) — Ftr (i))/m;

aa (i)=(ax (i)^2 +ay (i)^2)^0.5;

for i=1:n

t (i+1)=t (i)+dt;

%Скорости

Vy (i+1)=Vy (i)+(P (i)*sin (Phi (i)) — m*g (i) — X (i)*sin (teta (i))+Y (i)*cos (teta (i)))*dt/m;

if b>=1

Vy (i+1)=0;

end;

Vx (i+1)=Vx (i)+(P (i)*cos (Phi (i)) — X (i)*cos (teta (i)) — Y (i)*sin (teta (i)) — Ftr (i))*dt/m;

V (i+1)=((Vx (i+1)^2)+Vy (i+1)^2)^0.5;

%Угловая скорость

Phi (i+1)=Phi (i);

teta (i+1)=atan (Vy (i+1)/Vx (i+1));

Alpha (i+1)=Phi (i+1) — teta (i+1);

%Координаты

y (i+1)=y (i)+Vy (i+1)*dt;

x (i+1)=x (i)+Vx (i+1)*dt;

%Ускорение свободного падения

g (i+1)=g0*6 300 000² /(6 300 000+y (i+1))^2;

%Плотность

po (i+1)=p0*(1-y (i+1)/44 300)^5.236;

%Аэродинамические коэффициенты и тяга

P (i+1)=PF*po (i+1)/p0;

Cy (i+1)=polyval (p1, Alpha (i+1)*180/pi)+0.6*polyval (p1, Alpha (i+1)*180/pi);

Cxi (i+1)=0.056 *Cy (i+1)^2;

Cx (i+1)=polyval (p2, Alpha (i+1)*180/pi)+ Cxi (i+1);

%Аэродинамические силы

X (i+1)=po (i+1)*(V (i+1)^2)*Cx (i+1)*S/2;

Y (i+1)=po (i+1)*(V (i+1)^2)*Cy (i+1)*S/2;

%Сила трения

N (i+1)=m*g (i+1) — Y (i+1)*cos (teta (i+1))+X (i+1)*sin (teta (i+1)) — P (i+1)*sin (Phi (i+1));

Ftr (i+1)=N (i+1)*f;

%Ускорения

ay (i+1)=(P (i+1)*sin (Phi (i+1)) — m*g (i+1) — X (i+1)*sin (teta (i+1))+Y (i+1)*cos (teta (i+1)))/m;

ax (i+1)=(P (i+1)*cos (Phi (i+1)) — X (i+1)*cos (teta (i+1)) — Y (i+1)*sin (teta (i+1)) — Ftr (i+1))/m;

aa (i+1)=(ax (i+1)^2 +ay (i+1)^2)^0.5;

if Ftr (i+1)<=0

b=0;

end;

if b<=0

Ftr (i+1)=0;

end;

PhiGrad (i+1)=Phi (i+1)*180/pi;

AlphaGrad (i+1)=Alpha (i+1)*180/pi;

end;

Приложение 2

%Начальные данные

m=4820; %масса самолёта

p0=1.25; %плотность при y=0

f=0.02; %коэффициент трения сухого грунта

dt=0.01; %шаг интегрирования

Iz=5500; %момент инерции

Lp=9.23; %расстояние между центром масс и центром давления руля

Sp=10; %площадь руля

Ss=2; %площадь стабилизатора

Ls=9.23; %расстояние между центром масс и центром давления стабилизатора

delr0=-20*pi/180;%угол отклонения рулей

g0=9.8; %ускорение свободного падения при y=0

S=69.2; %площадь крыла

V0=0; %начальная скорость

n=10 000; %количество итераций

i=1; %счётчик

F (i)=0;

k0 (i)=0;

k1 (i)=0;

k2 (i)=0;

t (i)=0; %начальное время

y (i)=0; %начальная высота

x (i)=0; %начальное положение самолёта

teta (i)=0; %угол между вектором скорости и и горизонтом

w (i)=0; %Угловая скорость

delr (i)=delr0; %Угол отклонения рулей

Phi (i)=8*pi/180;%Угол тангажа

AlphaZV (i)=0;

b=1;

%Решение

%Ускорение свободного падения

g (i)=g0*6 300 000² /(6 300 000+y (i))^2;

%Проекции начальной скорости на оси OX и OY

Vx (i)=V0*cos (teta (i));

Vy (i)=V0*sin (teta (i));

V (i)=((Vx (i)^2)+Vy (i)^2)^0.5;

%Углы

Alpha (i)=Phi (i) — teta (i);

AlphR (i)=delr (i) — Alpha (i);

%Плотность

po (i)=p0*(1-y (i)/44 300)^5.236;

%Тяга

P (i)=0.9*m*g (i)*po (i)/p0;

%Скоростной напор

q (i)=0.5*po (i)*V (i)^2;

%Аппроксимация аэродинамических коэффициентов

cyPoints=[-0.84 -1.02 -0.75 -0.7 -0.63 -0.4 0.06 0.4 0.63 0.7 0.75 1.02 0.84 0.5];

cxPoints=[-0.2 -0.11 -0.08 -0.05 -0.04 -0.04 0.015 0.02 0.04 0.05 0.8 0.11 0.2 0.5];

cyRulPoints=[-0.78 -0.97 -0.78 -0.64 -0.57 -0.34 0 0.34 0.57 0.64 0.69 0.97 0.78 0.44];

AlphaPoints=[-20 -16 -13 -10 -8 -5 0 5 8 10 13 16 20 25];

cyzvPoints=[-0.78 -0.97 -0.78 -0.64 -0.57 -0.34 0 0.34 0.57 0.64 0.69 0.97 0.78 0.44];

p1=polyfit (AlphaPoints, cyPoints, 3);

p2=polyfit (AlphaPoints, cxPoints, 2);

p3=polyfit (AlphaPoints, cyRulPoints, 3);

p4=polyfit (AlphaPoints, cyzvPoints, 3);

%Аэродинамические коэффициенты

Cy (i)=1.5*polyval (p1, Alpha (i)*180/pi);

Cxi (i)=0.056*Cy (i)^2;

Cx (i)=1.18*polyval (p2, Alpha (i)*180/pi)+ Cxi (i);

Cn (i)=Cy (i)*cos (Alpha (i))+Cx (i)*sin (Alpha (i));

Cyzv (i)=polyval (p4, AlphaZV (i)*180/pi);

Cxizv (i)=0.0626*Cyzv (i)^2;

Cxzv (i)=polyval (p2, AlphaZV (i)*180/pi)+Cxizv (i);

Cnzv (i)=Cyzv (i)*cos (AlphaZV (i))+Cxzv (i)*sin (AlphaZV (i));

CyR (i)=polyval (p3, AlphR (i)*180/pi);

CxiR (i)=0.0626*CyR (i)^2;

CxR (i)=polyval (p2, AlphR (i)*180/pi)+CxiR (i);

CnR (i)=CyR (i)*cos (AlphR (i))+CxR (i)*sin (AlphR (i));

Cys (i)=polyval (p3, Alpha (i)*180/pi);

Cnss (i)=Cys (i)*cos (Alpha (i))+Cx (i)*sin (Alpha (i));

%Аэродинамические силы

X (i)=q (i)*Cx (i)*S;

Y (i)=q (i)*Cy (i)*S;

%Сила реакции опоры

N (i)=m*g (i) — Y (i)*cos (teta (i))+X*sin (teta (i)) — P (i)*sin (Phi (i));

%Сила трения

Ftr (i)=(N (i))*f;

PhiGrad (i)=Phi (i)*180/pi;

AlphaGrad (i)=Alpha (i)*180/pi;

AlphRGrad (i)=AlphR (i)*180/pi;

delrGrad (i)=delr (i)*180/pi;

AlphaZV1 (i)=0;

Cyzv1 (i)=0;

Cxizv1 (i)=0;

Cxzv1 (i)=0;

AlphaZV2 (i)=0;

Cyzv2 (i)=0;

Cxizv2 (i)=0;

Cxzv2 (i)=0;

Cnzv1 (i)=0;

Cnzv2 (i)=0;

%Аэродинамические моменты

Mss (i)=q (i)*Lp*(Ss*Cnss (i));

Mds (i)=q (i)*Lp*Ss*Cnzv (i);

MsR (i)=q (i)*Lp*Sp*CnR (i);

for i=1:n

t (i+1)=t (i)+dt;

Vy (i+1)=Vy (i)+(P (i)*sin (Phi (i)) — m*g (i) — X (i)*sin (teta (i))+Y (i)*cos (teta (i)))*dt/m;

if b>=1

Vy (i+1)=0;

end;

Vx (i+1)=Vx (i)+(P (i)*cos (Phi (i)) — X (i)*cos (teta (i)) — Y (i)*sin (teta (i)) — Ftr (i))*dt/m;

V (i+1)=((Vx (i+1)^2)+Vy (i+1)^2)^0.5;

%Угловая скорость

if b<=0

F (i)=q (i)*Lp*(-Ss*Cnss (i)+Ss*Cnzv (i)+Sp*CnR (i))/Iz;

k0 (i)=dt*F (i);

AlphaZV1 (i)=atan (Lp*(w (i)+0.5*k0 (i))/V (i));

Cyzv1 (i)=polyval (p4, AlphaZV1 (i)*180/pi);

Cxizv1 (i)=0.0626*Cyzv1 (i)^2;

Cxzv1 (i)=polyval (p2, AlphaZV1 (i)*180/pi)+Cxizv1 (i);

Cnzv1 (i)=Cyzv1 (i)*cos (AlphaZV1 (i))+Cxzv1 (i)*sin (AlphaZV1 (i));

k1 (i)=dt*q (i)*Lp*(-Ss*Cnss (i)+Ss*Cnzv1 (i)+Sp*CnR (i))/Iz;

AlphaZV2 (i)=atan (Lp*(w (i) — k0 (i)+2*k1 (i))/(V (i)));

Cyzv2 (i)=polyval (p4, AlphaZV2 (i)*180/pi);

Cxizv2 (i)=0.0626*Cyzv2 (i)^2;

Cxzv2 (i)=polyval (p2, AlphaZV2 (i)*180/pi)+Cxizv2 (i);

Cnzv2 (i)=Cyzv2 (i)*cos (AlphaZV2 (i))+Cxzv2 (i)*sin (AlphaZV2 (i));

k2 (i)=dt*q (i)*Lp*(-Ss*Cnss (i)+Ss*Cnzv2 (i)+Sp*CnR (i))/Iz;

w (i+1)=w (i)+(k0 (i)+4*k1 (i)+k2 (i))*0.166;

else

F (i)=0;

k0 (i)=0;

AlphaZV1 (i)=0;

Cyzv1 (i)=0;

Cxizv1 (i)=0;

Cxzv1 (i)=0;

Cnzv1 (i)=0;

k1 (i)=0;

AlphaZV2 (i)=0;

Cyzv2 (i)=0;

Cxizv2 (i)=0;

Cxzv2 (i)=0;

Cnzv2 (i)=0;

k2 (i)=0;

w (i)=0;

w (i+1)=0;

end;

Phi (i+1)=Phi (i)+w (i)*dt; %угол тангажа

teta (i+1)=atan (Vy (i+1)/Vx (i+1));%угол между вектором скорости и горизонтом

delr (i+1)=delr0; %угол отклонения руля высоты

Alpha (i+1)=Phi (i+1) — teta (i+1); %угол атаки

AlphR (i+1)=delr (i+1)+Phi (i+1) — teta (i+1); %угол атаки руля высоты

AlphaZV (i+1)=atan (w (i+1)*Lp/V (i+1)); %приращение угла атаки стабилизатора

%Координаты

y (i+1)=y (i)+Vy (i+1)*dt;

x (i+1)=x (i)+Vx (i+1)*dt;

%Ускорение свободного падения

g (i+1)=g0*6 300 000² /(6 300 000+y (i+1))^2;

%Плотность и тяга

po (i+1)=p0*(1-y (i+1)/44 300)^5.236;

P (i+1)=0.9*m*g (i)*po (i+1)/p0;

%Аэродинамические коэффициенты

Cy (i+1)=1.5*polyval (p1, Alpha (i+1)*180/pi);

Cxi (i+1)=0.056*Cy (i+1)^2;

Cx (i+1)=1.18*polyval (p2, Alpha (i+1)*180/pi)+ Cxi (i+1);

Cn (i+1)=Cy (i+1)*cos (Alpha (i+1))+Cx (i+1)*sin (Alpha (i+1));

CyR (i+1)=polyval (p3, AlphR (i+1)*180/pi);

CxiR (i+1)=0.0626*CyR (i+1)^2;

CxR (i+1)=polyval (p2, AlphR (i+1)*180/pi)+CxiR (i+1);

CnR (i+1)=CyR (i+1)*cos (AlphR (i+1))+CxR (i+1)*sin (AlphR (i+1));

Cyzv (i+1)=polyval (p4, AlphaZV (i+1)*180/pi);

Cxizv (i+1)=0.0626*Cyzv (i+1)^2;

Cxzv (i+1)=polyval (p2, AlphaZV (i+1)*180/pi)+Cxizv (i+1);

Cnzv (i+1)=Cyzv (i+1)*cos (AlphaZV (i+1))+Cxzv (i+1)*sin (AlphaZV (i+1));

Cys (i+1)=polyval (p3, Alpha (i+1)*180/pi);

Cnss (i+1)=Cys (i+1)*cos (Alpha (i+1))+Cx (i+1)*sin (Alpha (i+1));

%Аэродинамические силы

q (i+1)=0.5*po (i+1)*V (i+1)^2;

X (i+1)=q (i+1)*Cx (i+1)*S;

Y (i+1)=q (i+1)*Cy (i+1)*S;

%Сила трения

N (i+1)=m*g (i+1) — Y (i+1)*cos (teta (i+1))+X (i+1)*sin (teta (i+1)) — (i+1)*sin (Phi (i+1));

Ftr (i+1)=N (i+1)*f;

if Ftr (i+1)<=0

b=0;

end;

if b<=0

Ftr (i+1)=0;

end;

PhiGrad (i+1)=Phi (i+1)*180/pi;

AlphaGrad (i+1)=Alpha (i+1)*180/pi;

AlphRGrad (i+1)=AlphR (i+1)*180/pi;

delrGrad (i+1)=delr (i+1)*180/pi;

Mss (i+1)=q (i+1)*Lp*Ss*Cnss (i+1);

Mds (i+1)=q (i+1)*Lp*Ss*Cnzv (i+1);

MsR (i+1)=q (i+1)*Lp*Sp*CnR (i+1);

end