Проблемы машиностроения и надежности машин, 2021, № 4, стр. 70-80

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

В то же время, в России и за рубежом проектируются летательные аппараты с полетными числами Маха М = 6…15, внешние обводы которых имеют сложную форму, и обтекание которых может сопровождаться указанными выше химическими превращениями. Для достижения больших скоростей полета конструкция таких аппаратов должна иметь высокое весовое совершенство, которое обеспечит размещение на гиперзвуковых летательных аппаратах (ГЗЛА) необходимого для решения боевых задач бортового оборудования и боевой части. Конструкции ГЗЛА могут иметь тонкостенные панели обшивки и упругое оперение, и, следовательно, они могут быть подвержены флаттеру и возможным разрушениям от него. Расчеты границ устойчивости к флаттеру таких аппаратов на основе упрощенных методик могут давать количественно неверные предсказания границ устойчивости [2]. Таким образом, высока актуальность разработки более совершенной методики расчета аэроупругой устойчивости, учитывающей как реальную геометрию внешних обводов тела, так и химические превращения в воздухе, и создания на ее основе программных модулей и библиотек для расчета флаттера в сверхзвуковом потоке.

Методика расчета аэроупругой устойчивости. Выведем систему уравнений движения аэроупругой системы в обобщенных координатах из системы дифференциальных уравнений движения Лагранжа [3]

где L – лагранжиан, равный разности кинетической и потенциальной энергий; $Q_{i}^{n}$ – непотенциальные обобщенные силы, вызванные аэродинамическим воздействием; ${{u}_{i}}$ – обобщенные координаты, в качестве которых примем амплитуды разложения вектора перемещений поверхности упругого тела по собственным формам колебаний (т.е., по сути, используется метод Бубнова–Галеркина). Будем считать, что известно статическое аэроупругое положение равновесия, которое далее исследуется на устойчивость.

Для получения выражений для линеаризованных обобщенных сил рассмотрим элементарную работу ${{\delta }}W$, совершаемую аэродинамическим давлением на элементарном приращении обобщенных координат ${{\delta }}{{u}_{j}}$

где ${{\vec {p}}_{n}}$ – приращение вектора давления, действующего на тело (вязкими напряжениями в воздухе будем пренебрегать); S – площадь поверхности тела; $~{{\vec {w}}_{i}}\left( {x,y,z} \right)$ – распределение перемещений по i-й собственной форме; ${{\vec {w}}_{j}}{{\delta }}{{u}_{j}}$ – элементарное движение тела при вариации обобщенной координаты; N – число учитываемых собственных форм.

Из вариационного принципа Даламбера–Лагранжа [3] следует, что

Таким образом, вектор-столбец линеаризованных обобщенных аэродинамических сил можно представить в виде произведения матрицы

называемой матрицей аэродинамической жесткости, и вектора-столбца обобщенных координат. Здесь p – возмущение давления, вызванное деформацией тела по j-й собственной форме; $\vec {n}$ – нормаль к поверхности тела, направленная внутрь тела. Поскольку в статье рассматриваются высокие сверхзвуковые и гиперзвуковые скорости и относительно небольшие объекты, то характерное время движения частицы газа вдоль тела на порядок меньше характерных периодов колебаний тела по собственным формам, которые имеет смысл учитывать в расчетах на флаттер (малые числа Струхаля). Следовательно, обтекание можно считать квазистационарным, и распределение давления зависит лишь от перемещения, но не от скорости движения поверхности тела. Матрицей аэродинамического демпфирования, таким образом, пренебрегается, что дает несколько меньшие критические скорости флаттера и идет “в запас”.

Выражения для кинетической T и потенциальной U энергии линейно упругого тела в обобщенных координатах имеют вид

что дает выражение для лагранжиана L = T – U. В этих выражениях матрица $\hat {M}$ – матрица масс; матрица $\hat {K}$ – матрица конструкционной жесткости. Уравнения Лагранжа движения тела в матричной форме примут вид

Здесь также учтена матрица конструкционного демпфирования $\hat {D}$, которую можно задать при известных значениях коэффициентов демпфирования по каждому собственному колебанию.

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

Уравнение неразрывности $\frac{{\partial {{\rho }}}}{{\partial t}} + \nabla \cdot \left( {{{\rho }}U} \right) = 0$.

Уравнения Навье–Стокса $\frac{{\partial (\rho U)}}{{\partial t}}\, + \,\nabla \, \cdot \,(\rho U\, \otimes U) = - \nabla p\, + \,\nabla \, \cdot \,\tau $, где ${{(U \otimes U)}_{{ij}}}$ = ${{(U{{U}^{T}})}_{{ij}}}$ = = ${{U}_{i}}{{U}_{j}}$, и компоненты тензора τ выражаются как ${{\tau }_{{ij}}}$ = ${{\mu }}\left( {{{\nabla }_{j}}{{U}_{i}} + {{\nabla }_{i}}{{U}_{j}} - \frac{2}{3}{{{{\delta }}}_{{ij}}}\left( {\nabla \cdot U} \right)} \right)$, а ${{{{\delta }}}_{{ij}}}$ – символ Кронекера.

Уравнение энергии, записанное через полную энтальпию

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

Получив уравнения движения аэроупругой системы и вычислив матрицу ${{\hat {K}}_{a}}$, решается задача на комплексные собственные значения. Представляя u как $u = C{{e}^{{{{\lambda }}t}}}$, находим

В результате получаем набор N собственных значений λ. Критерием устойчивости является условие $\operatorname{Re} ~\lambda \leqslant 0$. При наличии вещественного положительного собственного значения имеется дивергенция, комплексного собственного значения с положительной вещественной частью – флаттер.

Практическая реализация методики. В настоящей статье описывается методика расчета флаттера, при которой обобщенные аэродинамические силы вычисляются на основе численного моделирования обтекания (CFD) реальной геометрии конструкции с возможностью учета происходящих в воздухе химических реакций. Как показано, в предыдущем разделе, в силу больших скоростей полета и относительно небольших размеров рассматриваемых объектов, предполагается, что обтекание колеблющейся конструкции квазистационарно. Конкретная реализация методики основана на использовании двух коммерческих программных пакетов, используемых в инженерной практике: MSC.Nastran [5, 6] – для расчета собственных частот и форм колебаний конструкции в пустоте и для решения задачи на собственные значения связанной аэроупругой системы и Ansys CFX [7] – для решения задачи о квазистационарном обтекании конструкции. Для интеграции этих двух программных пакетов, вычисления обобщенных аэродинамических сил, формирования матрицы аэродинамической жесткости и обработки результатов разработаны дополнительные программные модули.

Блок-схема проведения расчетов показана на рис. 1. Собственные формы и частоты колебаний (как корпусные, так и рулевых поверхностей) рассчитываются в MSC.Nastran стандартным методом конечных элементов (решатель SOL 103) и с использованием разработанного программного модуля передаются в газодинамический пакет Ansys CFX. В нем рассчитывается обтекание конструкции при деформациях корпуса по собственным формам. Далее, используя второй разработанный модуль, вычисляются обобщенные аэродинамические силы, формирующие матрицу аэродинамической жесткости. Эта матрица передается в MSC.Nastran для расчета комплексных собственных частот колебаний связанной аэроупругой системы (решатель SOL 110), который проводится методом Бубнова–Галеркина (в терминологии Nastran – “в модальных координатах”), что требует сравнительно небольшого числа учитываемых собственных форм при вычислении обобщенных аэродинамических сил. Согласно теории, критерием неустойчивости (дивергенции или флаттера) является наличие собственной частоты с положительной вещественной частью.

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

На рис. 2а, б представлена геометрия первого исследуемого объекта – беспилотного маневренного ЛА. Все расчеты для данной модели были проведены для параметров полета на высоте 25 км и числе Маха М = 5. Описание низших собственных форм колебаний модели приведено в табл. 1 и на рис. 3. В расчетах учитывались две корпусные собственные формы колебаний (первая и вторая изгибные формы с номерами 8 и 10), а также трех форм колебаний рулей (изгибная, изгибно-крутильная и вращательная – номера 20, 28, 30). Стоит отметить, что в табл. 1 несколько частот соответствуют одной форме, поскольку соответствующая собственная форма реализуется на различных рулях, а малое отличие в значениях частот обусловлено несимметричным расположением бортового оборудования и полезной нагрузки.

Расчет комплексных собственных значений в модуле SOL110 MSC.Nastran показал отсутствие флаттера (табл. 2).

Об этом свидетельствуют нулевые значения аэродинамического демпфирования $\Delta = - \operatorname{Re} \left( \lambda \right)$, приведенные для всех рассчитанных физических частот колебаний $\omega = \frac{{\operatorname{Im} \left( \lambda \right)}}{{2\pi }}$, которые практически не отличаются от частот колебаний в пустоте. Это вызвано достаточно большими собственными частотами колебаний рулей и большой высотой полета на рассмотренном режиме (25 км), при которой аэродинамическое воздействие потока на ракету мало. Результаты расчета подтверждены натурными испытаниями.

На рис. 2в, г представлены геометрии второго и третьего исследуемого объекта – оперение маневренного ЛА. Согласно экспериментальным данным, использование трапециевидного руля (рис. 2г) не вызывало флаттер на различных расчетных режимах полета. В то же время замена трапециевидной геометрии оперения на прямоугольную (рис. 2в) приводила к флаттеру. Расчеты по предлагаемой в этой статье методике подтверждают экспериментальные результаты. Все расчеты проводились с параметрами, соответствующими нулевой высоте над уровнем моря.

На рис. 4 изображены, а в табл. 3 описаны полученные собственные формы колебания второй модели в пустоте.

В табл. 4 приведены полученные по результатам расчета согласно методике результаты для оперения прямоугольной формы. Сначала был проведен расчет для числа Маха М = 2.5, по результатам которого был сделан вывод о флаттере связанного типа по второй форме на этом режиме полета. Далее были проведены расчеты для М = 2 и М = 1.5 (отметим, что в последнем случае предположение о квазистационарности потока заведомо неверно, и этот расчет представляет чисто академический интерес). Полученные значения коэффициентов аэродинамического демпфирования также свидетельствуют о наличии того же типа флаттера.

Аналогичные расчеты были проведены для трапециевидной модели оперения. На рис. 5 изображены, а в табл. 5 описаны полученные собственные формы модели. Сначала был проведен расчет для М = 2, затем на М = 6 (табл. 6). При этом Δ остался равен нулю, что указывает на отсутствие флаттера.

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

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