Известия РАН. Теория и системы управления, 2021, № 6, стр. 24-34

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

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

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

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

Второй подход предложен в [5, 6] и основан на знании координаты точки отрыва потока и разрушении вихревой системы на аэродинамическом профиле, которая в летном эксперименте не наблюдается. Применительно к нестационарным аэродинамическим характеристикам самолета этот подход рассмотрен в [7, 8]. По существу в этом случае самолет отождествляется с некоторым аэродинамическим профилем с эквивалентными аэродинамическими характеристиками.

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

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

1. Математическая модель нестационарных аэродинамических характеристик самолета. Рассмотрим математическую модель только коэффициента подъемной силы летательного аппарата в целом в условиях нестационарного обтекания. Пусть выполнены следующие предположения:

• аэродинамические характеристики самолета в установившемся движении известны с высокой точностью;

• на интервале обработки полетных данных конфигурация самолета остается постоянной;

• аппроксимация нестационарных аэродинамических коэффициентов зависит от постоянных параметров, значения которых можно уточнить по полетным данным;

• для описания нестационарных аэродинамических коэффициентов можно применить метод аэродинамических переходных функций [3].

Пусть зависимость коэффициента подъемной силы самолета ${{c}_{y}}$ от параметров полета (угла атаки α, числа M, угловой скорости тангажа ${{\omega }_{z}}$) и параметров конфигурации на неустановившемся режиме полета может быть представлена в виде

где ${{c}_{{y\infty }}}(\alpha ,M,{{\delta }_{{{\text{конф}}}}})$ – коэффициент подъемной силы, полученный в аэродинамических продувках или расчетным путем с учетом полетной конфигурации самолета; $\Delta {{c}_{y}}(t,\alpha ,M,{{\delta }_{{{\text{конф}}}}};\dot {\alpha },{{\dot {\omega }}_{z}})$ – поправка коэффициента подъемной силы в условиях неустановившегося движения с учетом полетной конфигурации самолета; ${{\delta }_{{{\text{конф}}}}} = ({{\delta }_{{\text{в}}}},{{\delta }_{{{\text{пр}}}}},{{\delta }_{{{\text{закр}}}}},{{\delta }_{{\text{и}}}},{{\delta }_{{{\text{т}}{\text{.щ}}}}},{{\delta }_{{\text{ш}}}},{\text{и}}\;{\text{т}}.{\text{д}}.)$ – вектор переменных, определяющих полетную конфигурацию самолета (отклонение руля высоты, предкрылков, закрылков, интерцепторов, тормозных щитков, шасси и т.д).

В силу предположения а) следует ожидать, что поправка $\Delta {{c}_{y}}(t,\alpha ,M,{{\delta }_{{{\text{конф}}}}};\dot {\alpha },{{\dot {\omega }}_{z}})$ мала по абсолютному значению и может быть представлена в виде

где $\Delta {{c}_{{y\infty }}}(\alpha ,M,{{\delta }_{{{\text{конф}}}}})$ – поправка коэффициента подъемной силы в установившемся движении; $\Delta {{c}_{{y\;{\text{дин}}}}}(t,\dot {\alpha },{{\dot {\omega }}_{z}})$ – динамическая поправка коэффициента подъемной силы, обусловленная нестационарными аэродинамическими эффектами в неустановившемся движении.

Поправку $\Delta {{c}_{{y\infty }}}(\alpha ,M,{{\delta }_{{{\text{конф}}}}})$ удобно представить в виде суммы некоторых функций от наиболее значимых переменных с неизвестными коэффициентами:

где ${{g}_{i}}(\alpha ,M,{{\delta }_{{{\text{конф}})}}}$ – известные функции и их количество N_g, которые подбираются в процессе решения задачи идентификации; ${{c}_{1}},{{c}_{2}}, \ldots $ – неизвестные коэффициенты, значения которых оцениваются по полетным данным на неустановившемся маневре. Примером суммы (1.1) может служить выражение $\Delta {{c}_{{y\infty }}}(\alpha ,M,{{\delta }_{{{\text{конф}}}}}) = {{c}_{1}}\alpha + {{c}_{2}}{{\delta }_{{\text{в}}}}.$

Для определения поправки $\Delta {{c}_{{y\;{\text{дин}}}}}(t,\dot {\alpha },{{\dot {\omega }}_{z}})$ применим метод аэродинамических переходных функций. Рассмотрим $\Delta {{c}_{{y\,{\text{дин}}}}}$ вида

где ${{h}_{\alpha }},{{h}_{{{{\omega }_{z}}}}}$ – аэродинамические переходные функции; b_a – средняя аэродинамическая хорда; $V$ – скорость полета.

Выражение (1.2) в уравнения динамики продольного движения самолета впервые введено в работе [2] и является результатом эвристических рассуждений о структуре аэродинамических коэффициентов как функционалов от угла атаки $\alpha $ и угловой скорости тангажа ${{\omega }_{z}}$. Однако ни из теории нестационарного обтекания, ни из аэродинамического эксперимента не ясно как задавать аэродинамические переходные функции в аналитическом виде на конкретном маневре самолета. Поэтому для выявления вида аэродинамических переходных функций и накопления статистической информации необходимо идентифицировать ${{h}_{\alpha }},{{h}_{{{{\omega }_{z}}}}}$ вместе с другими аэродинамическими поправками на различных тестовых маневрах в натурном эксперименте. В такой постановке задача идентификации нестационарных аэродинамических характеристик остается все еще сложной. Однако, учитывая достаточно общий характер структуры уравнения (1.2), заменим его выражением

Интегралы в (1.3) имеют одинаковую структуру и могут быть изучены как один интеграл с переменным верхним пределом

где $u(\tau ) = (\alpha (\tau )\;{\text{или}}\;{{\omega }_{z}}(\tau ))$.

Продифференцируем $\Delta C(t)$ по t, получим

Выберем переходную функцию такой, что

где $a = {\text{const}}$. Тогда совокупность выражений (1.4)–(1.6) может быть записана в виде обыкновенного дифференциального уравнения первого порядка относительно функции $\Delta C(t)$: с нулевым начальным условием $\Delta C(0) = 0$. В этом уравнении выполнена замена обозначений $h(0;u(t)) = \varphi (u(t))$. Предположим, что на интервале наблюдений. По условиям устойчивости решений уравнения (1.7) рассматриваются только отрицательные значения параметра a < 0.

Решая уравнение (1.7) с нулевым начальным условием и входным сигналом $\phi (u(\tau ))\dot {u}(\tau )$, получим

Функция $\varphi (u(\tau ))$ является неизвестной. Предположим, что на интервале $[0,T]$ функция $\varphi (u(t))$ допускает аппроксимацию ее полиномом по $u(t)$ с неизвестными коэффициентами:

где $K_{m}^{{(u)}}$ – неизвестные коэффициенты полинома; M – наивысшая степень полинома, задается исследователем в процессе решения задачи идентификации.

Подстановка функции $\varphi (u(t))$ в (1.8) приводит к выражению

Вычислим интеграл в (1.9) по частям, получим

где что является решением линейного дифференциального уравнения

В отличие от (1.9) выражение (1.10) не содержит производной процесса $u(t)$, которая либо не наблюдается в натурном эксперименте, либо вычисляется с большими погрешностями. Поэтому выражение (1.10) является предпочтительным для получения различных расчетных соотношений.

Выполняя в (1.10) замену переменной $u(t)$ на $\alpha (t)$ и ${{\omega }_{z}}(t)$, получим систему выражений, описывающих динамические поправки (1.2) коэффициента аэродинамической подъемной силы:

где все обозначения понятны из предыдущего текста. Здесь параметры ${{a}_{\alpha }},{{a}_{\omega }},K_{m}^{{(\alpha )}},K_{m}^{{(\omega )}},{{M}_{\alpha }}$ и M_ω неизвестны и должны определяться в процессе решения задачи идентификации по полетным данным.

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

Для определения неизвестных параметров в выражениях (1.1) и (1.11) по полетным данным применим частотно-временной метод [8] (см. Приложение). Вычислим на множестве частот $\Omega $ финитные преобразования Фурье F_T (с учетом обозначений, указанных ниже) всех составляющих поправок коэффициента аэродинамической подъемной силы:

где $\Delta {{c}_{{y1}}}(T),\Delta {{c}_{{y2}}}(T)$ – значения $\Delta {{c}_{{y1}}}(t),\Delta {{c}_{{y2}}}(t)$ в конечный момент времени T.

Таким образом, вектор неизвестных параметров $\theta $ имеет вид

Множество допустимых значений параметров определяется условиями устойчивости решений дифференциальных уравнений ${{a}_{\alpha }} < 0,\;{\kern 1pt} {{a}_{\omega }} < 0$, остальные параметры могут принимать любые значения. В зависимости от полетной конфигурации некоторые параметры могут отсутствовать.

Сформируем невязку и критерий метода наименьших квадратов (МНК) в частотной области

где ${{\lambda }_{\alpha }},{{\lambda }_{\omega }},{{\lambda }_{i}}$ – весовые коэффициенты, подбираются опытным путем; заранее верхний предел суммирования указать не представляется возможным, суммирование ведется по всем доступным частотам ω_k. Первые три слагаемых в критерии $J(\theta )$ служат для регуляризации решения задачи идентификации. Решение о включении регуляризирующих составляющих принимается в процессе решения задачи идентификации.

Задача идентификации неизвестных параметров сводится к задаче минимизации критерия МНК:

которая решается методами нелинейного математического программирования (например, программой fmincon математического пакета MATLAB).

На множестве частот $\Omega $ финитное преобразование Фурье постоянных, которые могут входить в составляющие коэффициента подъемной силы в установившемся движении самолета, равны нулю (см. ПРИЛОЖЕНИЕ). Поэтому для полной оценки всех постоянных дополнительно необходимо решить задачу идентификации во временной области на уже вычисленных оценках поправок $\Delta {{\hat {c}}_{{y\infty }}}(\alpha ,M,{{\delta }_{{{\text{конф}}}}})$ и $\Delta {{\hat {c}}_{{y\,{\text{дин}}}}}(t,\dot {\alpha },{{\dot {\omega }}_{z}})$.

2. Идентификация коэффициента подъемной силы ближнемагистрального самолета. Рассмотрим идентификацию коэффициента подъемной силы ближнемагистрального самолета по полетным данным. Тестовые полеты проводились по программе исследований устойчивости и управляемости самолета в процессе его выхода на большие углы атаки в различных полетных конфигурациях. Тестовый маневр заключался в выводе самолета на большие углы атаки при сохранении балансировки по моменту тангажа, т.е. ${{\dot {\omega }}_{z}} \sim 0$. Балансировка нарушалась только на коротких временных интервалах, когда самолет возвращался к условиям полета на малых углах атаки. Все тестовые маневры выполнялись на высотах примерно H = 5000 м и числах $M = 0.3 \cdots 0.4$. На рис. 1 показан пример полетных данных во взлетной конфигурации ${{\delta }_{{{\text{пр}}}}}$ = 24°; ${{\delta }_{{{\text{закр}}}}}$ = 16°. Далее рассмотрим определение математической модели только динамической поправки $\Delta {{c}_{{y\,{\text{дин}}}}}(t,\dot {\alpha },{{\dot {\omega }}_{z}})$. Измерялись все переменные движения и параметры полетной конфигурации. Полетный коэффициент подъемной силы $c_{y}^{ * }(\alpha ,M,{{\delta }_{{{\text{конф}}}}})$ (без учета тяги двигателя) вычислялся по формуле

где ${{c}_{y}}$ – коэффициент аэродинамической подъемной силы; ${{c}_{P}}$ – коэффициент тяги двигателей; φ_P – угол установки двигателей в плоскости $OXY$ связанной системы координат. Для каждой полетной конфигурации известны продувки в аэродинамической трубе (АДТ) и, следовательно, априорные значения коэффициента подъемной силы в установившемся движении ${{c}_{{y\infty }}}(\alpha ,M,{{\delta }_{{{\text{конф}}}}})$, которые можно принять за начальное приближение.

Требуется по измеренным данным, полученным в тестовом полете, уточнить значение коэффициента подъемной силы в установившемся движении и идентифицировать составляющую $\Delta {{c}_{{y\,{\text{дин}}}}}(t,\dot {\alpha },{{\dot {\omega }}_{z}})$, обусловленную нестационарными аэродинамическими эффектами.

Рассматривались следующие полетные конфигурации самолета: крейсерская, взлетная, заход на посадку, посадка. Поправку $\Delta {{c}_{{y\infty }}}(\alpha ,M,{{\delta }_{{{\text{конф}}}}})$ будем искать в виде

В процессе решения задачи идентификации были выбраны следующие порядки полиномов ${{M}_{\alpha }} = {{M}_{\omega }} = 4$; определялись оценки неизвестных параметров $({{c}_{1}}, \ldots ,{{c}_{4}},{{a}_{\alpha }},{{a}_{\omega }})$, коэффициентов полиномов $(K_{0}^{{(\alpha )}}, \ldots ,K_{4}^{{(\alpha )}};K_{0}^{{(\omega )}}, \ldots ,K_{4}^{{(\omega )}})$ и дополнительных параметров $(\Delta {{c}_{{y1}}}(T),\Delta {{c}_{{y2}}}(T))$; все углы и угловые скорости имеют размерности рад и рад/с соответственно.

Результаты идентификации коэффициента подъемной силы самолета по полетным данным с учетом конфигурации, статических и динамических поправок показаны на рис. 2–6. На рис. 2, а–6, а помечены следующие данные: полетные данные – сплошная линия светлого оттенка, вычисленная по формуле (2.1); оценка коэффициента подъемной силы в установившемся движении с учетом поправки $\Delta {{\hat {c}}_{{y0}}} + \Delta {{\hat {c}}_{{y\infty }}}(\alpha ,M,{{\delta }_{{{\text{конф}}}}})$ – пунктирная линия; оценка коэффициента подъемной силы $\hat {c}_{y}^{*}(\alpha ,M,{{\delta }_{{{\text{конф}}}}})$ в неустановившемся движении – сплошная линия. На рис. 2, б–6, б показаны оценки динамических поправок: сплошная линия – динамическая поправка $\Delta {{\hat {c}}_{{y1}}}(t,\dot {\alpha })$; пунктирная линия – динамическая поправка $\Delta {{\hat {c}}_{{y2}}}(t,{{\dot {\omega }}_{z}})$.

На рис. 2, а представлен коэффициент подъемной силы и его оценки в области малых углов атаки при безотрывном обтекании, полностью совпадающие с аэродинамическими продувками. На рис. 2, б даны оценки динамических поправок. Видно, что динамическая поправка $\Delta {{\hat {c}}_{{y2}}}(t,{{\dot {\omega }}_{z}})$ = 0, а поправка $\Delta {{\hat {c}}_{{y1}}}(t,\dot {\alpha })$ имеет порядок 10^–3 и ее влиянием на оценку коэффициента подъемной силы можно пренебречь.

На рис. 3, а–6, а показаны коэффициент подъемной силы и его оценки в диапазоне изменения углов атаки $0 < \alpha \leqslant 35^\circ $ с учетом полетной конфигурации. Разность сплошной и штриховой линий на рис. 3, а–6, а указывает на влияние нестационарных аэродинамических эффектов на аэродинамические характеристики самолета. Возврат самолета на малые углы атаки происходит с потерей подъемной силы. Наблюдается аэродинамический гистерезис (рис. 3, а; 4, а; 6, а). На рис. 5, а показана аэродинамическая нестационарность в виде автоколебаний по коэффициенту подъемной силы при гармоническом изменении угла атаки: автоколебания возникают как на малых, так и на больших углах атаки.

Соизмеримые динамические поправки, обусловленные $\dot {\alpha }$ и ${{\dot {\omega }}_{z}}$ приведены на рис. 3, б. На остальных рассмотренных тестовых режимах полета (см. рис. 4, б–6, б) динамические поправки, обусловленные ${{\dot {\omega }}_{z}}$, практически отсутствуют (т.е. на тестовом режиме самолет сбалансирован в продольном движении). Поэтому динамические поправки оценки коэффициента подъемной силы целиком определяются угловой скоростью угла атаки $\dot {\alpha }$.

В табл. 1 и 2 приведены оценки коэффициентов полиномов по $\alpha (t)$ и ${{\omega }_{z}}(t)$, упорядоченные по возрастанию степени полинома, и коэффициенты поправок, полученные на выбранных тестовых маневрах. Регуляризирущая добавка в критерии МНК имела вид ${{(0.1{{a}_{\alpha }})}^{2}} + {{(0.05{{a}_{\omega }})}^{2}}$.

Коэффициенты первой строки табл. 1 соответствуют динамической поправке $\Delta {{c}_{{y\,{\text{дин}}}}}(t,\dot {\alpha },\mathop {\dot {\omega }}\nolimits_z )$, которая является откликом линейной динамической системы только на входные сигналы $\dot {\alpha }(t)$ и ${{\dot {\omega }}_{z}}(t)$ [2, 3]. Это частный случай описания нестационарных аэродинамических эффектов при безотрывном обтекании [6], соответсвующее выражение имеет вид

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

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