Модель динамики полета одного из классов летательных аппаратов с учетом упругости конструкции

Civil Avition High TECHNOLOGIES

Vol. 19, № 01, 2016

УДК 629.7.015.4

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

А.И. НИКИТИН

Статья представлена доктором технических наук, профессором Вышинским В.В.

Остается актуальной проблема демпфирования возникающих упругих колебаний конструкции летательного аппарата (ЛА) средствами бортовых систем автоматического управления. При решении подобной задачи модель упругого ЛА является основой синтеза законов управления и анализа замкнутой системы "объект управления -регулятор". В общем случае задача математического моделирования динамики полета упругого ЛА распадается еще как минимум на две взаимосвязанные задачи, одна из которых - непосредственно моделирование поведения упругого ЛА, определяемого действующими на него силами, а другая - учет изменения действующих аэрогидродинамических сил и моментов, вызванного деформацией упругого ЛА и работой системы управления.

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

Ключевые слова: математическая модель, динамика полета, упругий летательный аппарат, система автоматического управления, ракета.

Исторически, первыми практическими задачами, решение которых потребовало учета упругости летательных аппаратов (ЛА) и его частей, были задачи определения динамических нагрузок и устойчивости конструкции, в том числе задачи о флаттере [1, 2]. Без учета деформации конструкции и действия инерционных сил решение подобных задач невозможно. По мере развития авиационной и ракетной техники происходило дальнейшее увеличение скоростей полета и относительное уменьшение жесткости конструкций ЛА. Стала актуальной проблема демпфирования (гашения) возникающих упругих колебаний конструкции средствами бортовых систем автоматического управления (САУ) [3, 4]. В этом случае модель упругого ЛА является основой синтеза законов управления и анализа замкнутой системы "объект управления - регулятор".

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

Наиболее точное решение первой задачи можно получить, используя метод конечных элементов (МКЭ). Однако такой подход имеет ряд недостатков. Во-первых, подготовка конечно-элементной модели - процесс достаточно трудоемкий. Во-вторых, решение задачи упругих деформаций с помощью МКЭ само по себе требует значительных вычислительных ресурсов, а расчет изменяющихся аэрогидродинамических сил, действующих на ЛА, (вторая задача) и необходимость учета отклонения управляющих поверхностей вызовут еще большее увеличение потребных ресурсов. В-третьих, при таком подходе невозможно решение задачи аналитического синтеза САУ гашения колебаний. По этим причинам наиболее рациональным подходом к решению части задач моделирования динамики упругого ЛА остаются упрощенные методы, основанные на замене реальной конструкции ЛА той или иной эквивалентной схемой.

В общем случае ЛА является упругим телом с распределенной массой и, согласно теории колебаний [1], имеет бесконечное число собственных форм (или тонов) и частот колебаний.

Vol. 19, № 01, 2016

Civil Avition High TECHNOLOGIES

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

Исходя из этого допущения, возможны две основные расчетные модели (или схемы): балочная (рис. 1,6) и в виде набора сосредоточенных масс (рис. 1,в). Первая хорошо подходит для описания движения по низшим тонам упругого ЛА с крылом большого удлинения. Число уравнений, составляющих математическую модель ЛА, в этом случае определяется числом учитываемых тонов (форм) колебаний и их типом (изгибные, крутильные). Вторая схема позволяет получить хороший эквивалент конструкции ЛА со стреловидным крылом малого удлинения, однако число уравнений, описывающих динамику конструкции ЛА, значительно возрастает и определяется числом используемых сосредоточенных масс.

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

Уравнение изгибных колебаний упругой балки под действием распределенной нагрузки с интенсивностью Р^,1;) (рис. 2) имеет вид [3]:

Civil Avition High TECHNOLOGIES

Vol. 19, № 01, 2016

где Е - модуль Юнга материала балки, I - момент инерции сечения балки относительно нейтральной оси, т - погонная масса балки.

Уравнение (1) допускает разделение переменных, поэтому его общее решение состоит из суммы частных и имеет вид:

оо

у в

1=1

В (2) Я. (2) - формы собственных колебаний, являющиеся решениями уравнения в частных производных

а Я^) - Т Н. нормальные координаты, подлежащие определению из системы обыкновенных дифференциальных уравнений (ОДУ) вида:

М Я : + М «.V = Р.

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

ь

М = I" Г2

и

Р = | Р^)^)^.

. 0

В исходном дифференциальном уравнении (1) изгибных колебаний балки не учтено рассеяние энергии за счет внутреннего трения в материале балки. В действительности это рассеяние всегда имеет место, благодаря чему собственные колебания затухают. В [3] указано, что рассеяние энергии можно учесть непосредственно в уравнении (4) слагаемым, пропорциональным Я.(1), а соответствующий коэффициент пропорциональности ^ (относительный коэффициент демпфирования) на практике обычно одинаков для всех тонов упругих колебаний. Принимая эту точку зрения, перепишем уравнения (4) реакции собственных форм на возмущающую силу в виде:

М.Я. + 2С.М.ю.я. + Мю.2я. = Р. С = 1,2,...,~).

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

Vol. 19, № 01, 2016

Civil Avition High TECHNOLOGIES

ных условий определяется конечное число собственных форм ^ (2) и частот изгибных колебаний балки по низшим тонам. Затем для этих форм вычисляются интегралы (5) и (6). Далее находится решение системы обыкновенных дифференциальных уравнений (7). После подстановки этого решения qi(t) и функций собственных форм в ряд (2), ограничиваясь, соответственно, рассмотрением конечного числа слагаемых, получим аппроксимацию искомого решения у(2, ^ .

Решения уравнения в частных производных (3) могут быть получены различными методами. На практике наибольшее распространение получил так называемый метод последовательных приближений, широко освещенный в литературе [2, 5]. Алгоритм метода был реализован в пакете программ МаШСаё. Тестирование расчетов проводилось путем сравнения полученных собственных форм и частот колебаний с известными теоретическими результатами для однородной балки постоянного сечения. Рассматривалась консольная и свободная неопертая балка. В обоих случаях рассчитанные и теоретические формы и частоты собственных колебаний совпали.

Полученное выше уравнение (7), описывающее изгибные колебания балки, фактически является математической моделью наиболее простого, с точки зрения схематизации, летательного аппарата - ракеты. С этого объекта и началась программная реализация модели динамики полета упругого ЛА. В [3] на основе анализа решения (7) показано, что нулевые тона колебаний неопертой балки описывают ее поведение в продольной плоскости, как твердого тела. Из общего уравнения (7) получены два частных: уравнение движения центра масс (ЦМ) в скоростной системе координат (СК) и уравнение вращения в продольной плоскости относительно ЦМ. Логично предположить, что модель упругой ракеты может быть получена дополнением традиционных уравнений динамики полета твердого ЛА в любой СК уравнениями колебаний балки по тонам выше нулевого.

Исходя из этого предположения, реализация компьютерной модели (в среде МАТЬАБ^тиНпк) заключалась в дополнении полной системы уравнений динамики полета ЛА в земной СК из [6] уравнениями изгибных колебаний (7). Для этого они были приведены к системе дифференциальных уравнений в форме Коши:

р 2  м     г

qJ=

где - скорости нормальных координат.

Проверку программной реализации, а также подтверждение выдвинутой ранее гипотезы, предполагалось осуществить на тестовом примере: используя исходные данные для ракеты, приведенные в [3], повторить моделирование и сравнить с имеющимися результатами.

Изменение действующих аэродинамических сил и моментов, вызванное деформацией упругого ЛА, в указанном примере осуществлялось упрощенно с использованием гипотезы стационарности. В этом случае аэродинамическая сила Уа (а значит и момент) однозначно

" a

*

определяется текущими значениями местных углов атаки а корпуса ракеты:

Ya = qDSj cya

0

где qD = 0.5рУ2 - скоростной напор, Б - площадь миделевого сечения, са(х) - погонная производная коэффициента подъемной силы по углу атаки как функция положения х элементарного участка корпуса ракеты ёх, Ь - длина ракеты.

Civil Avition High TECHNOLOGIES

Vol. 19, № 01, 2016

Остановимся более подробно на определении местных углов атаки. Рассмотрим некоторую точку О на осевой линии ракеты, находящуюся на расстоянии А! от ЦМ (рис. 3). Для удобства расположим в этой точке центр связанной с недеформированной ракетой системы координат хоу. Ракета движется относительно неподвижного воздушного потока со скоростью V. Угол атаки в точке О для ракеты, как для твердого невращающегося тела, определяется соот-

Г V ^

ношением а = arctg

V.

где Vx и Vy проекции вектора воздушной скорости на соответ-

x У

ствующие оси связанной СК. Для случая малых углов атаки можно принять:

а

- Vy/Vx - Vy/V

Вращение ракеты относительно ЦМ с угловой скоростью Ю2 приведет к появлению в точке О линейной скорости, проекция которой на ось оу будет иметь величину Ю2 А1 и направление (при положительном значении Ю2), противоположное направлению Vy . Также противоположное V направление будет иметь и скорость деформации (прогиба) ракеты в точке О, обозначенная у.

By

Осевая линия

после деформации ^

у

* "'У/

¡);М +.Г

а, "+

Рис. 3

Таким образом, с учетом вращения и скорости деформации ракеты изменится величина

проекции вектора скорости относительно воздушного потока на ось оу, а сам вектор, обозна-

*

ченный на рис. 3 как V , изменит свою величину и направление. Еще одним фактором, определяющим значение местного угла атаки в рассматриваемой точке О, является угол прогиба осе-

в этой точке. С учетом

„ „ Эу

вой линии ракеты, величина которого определяется производной —

Эх

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

* (Vy +wz Al + y) dy

а «----+ —

V Эх

wz Al y dy

а—z--— + —

V V Эх

Обобщая выражение (11) на все точки ракеты можно записать:

а*

V

V

Эх

О

О

Vol. 19, № 01, 2016

Civil Avition High TECHNOLOGIES

Кроме учета изменения аэродинамических сил от деформации в рассматриваемом примере учитывается изменение проекции вектора тяги Ут, вызванное дополнительным поворотом сопла двигателя вслед за осевой линией ракеты [3]:

(

Y = T

5 +

Эу Эх

Л

в у

где Т - тяга двигателя (рис. 3), 8 - угол поворота сопла,

Эу Эх

- угол прогиба осевой линии ра-

кеты в точке подвеса двигателя B. Выражения для аэродинамического момента и момента силы тяги получены умножением на соответствующие плечи сил (9) и (12).

Таким образом, сформирована математическая модель продольного движения упругой ракеты, учитывающая конечное число тонов изгибных колебаний и влияние этих колебаний на действующие силы и моменты. После подстановки числовых данных и функциональных зависимостей в эту модель и с учетом первых двух тонов изгибных колебаний в [3] получена упрощенная модель с постоянными коэффициентами в матричной форме для конкретных условий полета: скорости и высоты. Однако повторить этот переход от полных уравнений динамики полета к упрощенным уравнениям не представляется возможным, так как в [3] указаны не все исходные данные. Отсутствуют значение скорости V и высоты H полета, значение сечения миделя ракеты S и зависимость погонной производной коэффициента подъемной силы по углу атаки c^(x) . По этой причине, к сожалению, невозможно и повторить результаты [3] на реализованной компьютерной модели. Чтобы восполнить указанный недостаток исходных данных, было принято допущение о равномерном распределении погонной производной коэффициента

подъемной силы по углу атаки cy(x) вдоль ракеты: c^(x) = C^/l = const. С учетом этого допущения и опираясь на имеющиеся в [3] числовые данные модели (в матричной форме) были получены следующие значения величины SC^ = 553.7 м2/рад, воздушной скорости V = 345.8 м/с, а сама упрощенная модель продольного движения приняла вид:

- 0.381 s2 + 0.0887s

- 9.12 0.861s

x

a (s) " - 81

fl(s) -1.04

kx(s) 19.7

k2(s)_ 19.3

5

В матричном уравнении (13) Ф - угол тангажа, б = - символьное обозначение производной по времени. Для удобства это уравнение было преобразовано к общепринятому виду х = Ах + Ви и реализовано в компьютерной модели в виде одного стандартного 81шиНпк-блока. Дополняя компьютерную 81шиНпк-модель моделью сервопривода поворота камеры (сопла) двигателя 8 = 108и/(Б +10) и законом управления для входного сигнала сервопривода

8„ вида:

(лпп 3.23 9.94V ( 14.53 Y. ( 1ЛЛ„ 11.47 158 |п 5,, = I 4.22 +-+-k +1 -4.345--__ |q2 +1 -100.3--+-|й

s + 2.92 s

s + 0.745

s + 0.0293 s +1.925

получим замкнутую систему "объект - регулятор". Аналогичным образом была дополнена и модель на основе полной системы уравнений пространственного движения. Фактически, мо-

в

Civil Avition High TECHNOLOGIES

Vol. 19, № 01, 2016

дель упругой ракеты была реализована в среде МЛТЬЛВ/81шиНпк двумя различными способами. Результаты моделирования переходных процессов для обоих способов реализации представлены на рис. 4 и рис. 5.

а 04

0 03 002

-0 01

-0 02 -0 03

...... _Упрощенная модель ракеты с постоянными коэффициентами      Полная модель упругой ракеты

—Д. i

Чу»

3

i 1 i

5 t. е

2 16 t

0 5 О

•0 5 -1 -15 -2

_Упрощенная модель ракеты с постоянными коэффициентами * * * *  Полная мадепь упругой ракеты

;

¡Iii

; : i 1

! i i

I

:

1 i I i

! ! i i

5

1. с

Рис. 4

О 025 О 02 0016 0 01 0 005 О

-0.005 0.01 -0015 -С 02

_Упрощенная модель ракеты с постоянными коэффициентами

......1......

j

! f ] *

i

:

i I г

£ В t.c

002 0 015

оо1

0 045

=г О

-0 005 -0 01 -0 015 -0 02

_Упрощенная модель ракеты с постоянными коэффициентами

I j * Полная модель упругой ракеты

;

ß i

I i

i '}

vi

1 i ' 1

5

1.С

Рис. 5

Графики представляют собой реакцию по тангажу (рис. 4,а), управляющему сигналу сервопривода (рис. 4,6) и нормальным координатам (рис. 5) на одно и тоже возмущающее воздействие: импульс сигнала 8и величиной 1.72° и длительностью 0.2 с. На приведенных

рисунках сравниваются результаты моделирования упрощенной и полной модели упругой ракеты. Упрощенная модель в виде системы дифференциальных уравнений с постоянными коэффициентами (13) получена на основе уравнений (7) при рассмотрении нулевой и еще первых двух форм изгибных колебаний. Полная модель получена путем дополнения уравнений динамики полета твердого ЛА в земной СК из [6] уравнениями (8) для первых двух форм изгибных колебаний. Хорошее согласование графиков при различных способах реализации модели динамики свидетельствует о справедливости предположения, что уравнения гибкого ЛА могут быть получены дополнением уравнений динамики твердого ЛА уравнениями, описывающими дополнительные степени свободы (балки, либо сосредоточенной массы).

Сравнение переходных процессов также являлось дополнительной проверкой правильности реализации компьютерной модели.

Vol. 19, № 01, 2016

Civil Avition High TECHNOLOGIES

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

ЛИТЕРАТУРА

1. Гудков А.И., Лешаков П.С. Внешние нагрузки и прочность летательных аппаратов. -М.: "Машиностроение", 1968. - 470 с.

2. Бисплингхофф Р.Л., Эшли X., Халфмен Р.Л. Аэроупругость. - М.: Изд. иностранной литературы, 1958. - 799 с.

3. Кузовков Н.Т. Системы стабилизации летательных аппаратов (баллистических и зенитных ракет). - М.: "Высш. школа", 1976. - 304 с.

4. Кашин Г.М., Федоренко Г.И. Автоматическое управление продольным движением упругого самолета. - М.: "Машиностроение", 1974. - 312 с.

5. Гроссман Е.П., Пановко Я.Г. Упругие колебания частей самолета. - Л.: Ленинградская военно-воздушная инженерная академия, 1947. - 260 с.

6. Остославский И.В., Стражева И.В. Динамика полета. Траектории летательных аппаратов. - М.: "Машиностроение", 1969. - 500 с.

FLIGHT DYNAMICS MODEL OF ONE CLASS OF AIRCRAFT WITH A VIEW OF ELASTIC CONSTRUCTION

Nikitin A.I.

It remains urgent problem of damping of elastic vibrations occurring aircraft structure means the automatic control systems on board. In solving this problem the aircraft elastic model is the basis for the synthesis of control laws and analysis of closed-loop system "control object - the regulator." In general, the problem of mathematical modeling of flight dynamics of the elastic aircraft breaks for at least another two objectives, one of which - direct simulation of the behavior of elastic aircraft defined interacting forces, and the other - the account of the changes operating aerogidrodynamic forces and moments caused by the deformation elastic aircraft and work control systems.

This paper discusses the theoretical basis of the approach to the solution of this problem, based on the replacement of the actual design of aircraft by equivalent circuit and its implementation for the missiles, the most simple in terms of schematic, class of aircraft. At the same time accounting for changes in aerodynamic forces and moments caused by the elastic deformation of the aircraft, it was performed by help of stationary hypothesis.

Key words: mathematical model, flight dynamic, elastic aircraft, automatic control system, rocket.

REFERENCES

1. Gudkov A.I., Leshakov P.S. Vneshnie nagruzki i prochnost letatelnykh apparatov [External loads and strength of aircraft], Moscow, Mashinostroenie. 1968. 470 p.

2. Bisplinghoff R.L., Ashley H., Halfman R.L. Aerouprugost [Aeroelasticity], Moscow, Izdatelstvo inostrannoy literatury. 1958. 799 p.

3. Kuzovkov N.T. Sistemy stabilizatsii letatelnykh apparatov

4. Kashin G.M., Fedorenko G.I. Avtomaticheskoe upravlenie prodolnym dvizheniem upru-gogo samoleta [Elastic-structure aircraft longitudinal motion automatic control], Moscow, Mashi-nostroenie. 1974. 312 p.

Civil Avition High TECHNOLOGIES

Vol. 19, № 01, 2016

5. Grossman E.P., Panovko J.G. Uprugie kolebania chastey samoleta [Aircraft elastic parts oscillations], Leningrad, Leningradskaya voenno-vozdushnaya inzhenernaya akademia. 1947. 260 p.

6. Ostoslavsky I.V., Strazheva I.V. Dinamika poleta. Traektorii letatelnykh apparatov [Flight dynamics. Flight paths of aircrafts], Moscow, Mashinostroenie. 1969. 500 p.

СВЕДЕНИЯ ОБ АВТОРЕ

Никитин Александр Игоревич, 1981 г.р., окончил Таганрогский радиотехнический университет (2004), кандидат технических наук, инженер-конструктор ПАО "Таганрогский авиационный научно-технический комплекс им. Г.М. Бериева", автор 18 научных работ, область научных интересов - синтез законов управления, математическое моделирование, динамика полета летательных аппаратов, электронный адрес: savr.nik@gmail.com.