Численное моделирование конструкции на ранних стадиях проектирования вертолета Текст научной статьи по специальности «Строительство и архитектура»

УЧЕНЫЕ ЗАПИСКИ ЦАГИ

Том ХЫУ 2013 № 2

УДК 629.735.45.064

ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ КОНСТРУКЦИИ НА РАННИХ СТАДИЯХ

ПРОЕКТИРОВАНИЯ ВЕРТОЛЕТА

Е. В. КАСУМОВ

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

Ключевые слова: математическая модель, численный эксперимент, кинематический анализ, силовой анализ, напряженно-деформированное состояние, флаттер, вертолет, проектирование.

ВВЕДЕНИЕ

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

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

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

1. УРАВНЕНИЕ ДВИЖЕНИЯ КОНСТРУКЦИИ С УЧЕТОМ СИЛ ДЕМПФИРОВАНИЯ И ИЗМЕНЯЕМОЙ ПО ВРЕМЕНИ ВНЕШНЕЙ СИСТЕМЫ СИЛ

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

Му + СУ + КУ = Р, (1)

где С — матрица демпфирования, размерность которой совпадает с размерностью матриц К и М; М — матрица масс системы; К — матрица жесткости системы; Р — сосредоточенные, объемные и поверхностные силы, действующие на механизм.

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

С = аМ + ЬК. (2)

КАСУМОВ Евгений Владимирович

кандидат технических наук, ведущий инженер КГТУ им. А. Н. Туполева

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

Однако расчетчик не всегда в достаточной степени обладает экспериментальными данными. Одним из способов решения системы уравнения (1) является описанный в работе [1] метод модальных усечений, который давно применяется для моделирования динамического отклика элементов конструкции на действие внешней нагрузки. Метод реализован в таких широко распространенных расчетных конечно-элементных комплексах, как ADAMS, NASTRAN, ANSYS и давно используется для расчетов механизмов с учетом деформации звеньев.

Основная суть решения системы уравнения (1) методом модального усечения состоит в том, что линейные деформации тела в локальной системе координат рассматриваются в виде линейной комбинации конечного числа векторов ф (форма моды):

N

' = ^фг§г,

(3)

i=1

где g — амплитуда моды, 7 — номер моды.

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

Mv + Kv = 0.

(4)

Для каждого варианта граничных условий определяется набор собственных частот с применением конечно-элементной модели. Это в заданном расчетчиком количестве последовательность собственных (начиная с низшей) частот. Практика показывает, что для достаточной точности решения системы уравнения (1) необходимо, как правило, не менее 10 первых мод для каждого варианта граничных условий. Выбранный набор мод и называется модальным усечением.

В геометрическом смысле с учетом выражения (3) сложная форма колебания тела рассматривается как линейная комбинация форм собственных значений (рис. 1).

/

:gn

i=1

Рис. 1. Линейная комбинация форм собственных значений

В результате, матрица масс системы (1) имеет блочный вид. Ее блоки выражаются в десяти инвариантах инерции как зависимость матрицы масс от характеристик модального усечения:

Ыа = ;

Mtr =-A

J 2 + J3 gj

B;

Mrr = B

J7 -

Mtm = AJ

J8 + J8 г

3.

gj- Ji

J Vb^J

B;

(5)

M ™ = BJ

J4 + J gj

M = J

Здесь А и В — матрицы преобразования системы координат, g — модальные координаты. Инварианты инерции вычислены для узлов конечно-элементной модели, основанной на информации о массе каждого узла (см. [1]).

Мг Мт '

М = к Мгг М гт

_Мт МТт Мтт _

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

Обобщенная матрица жесткости К содержит блоки Ктт обобщенной матрицы жесткости структурных компонентов относительно модальных координат g, наряду с блоками тел, представленных абсолютно жесткими.

Поэтому обобщенная матрица жесткости имеет вид:

' Кп КГг Кт

К= К1 Кгг К гт

_ К1т К гт К тт

0 0 0 0 0 0 0 0 К

(7)

На основе сформированных матриц М, С, К, Р формируется основная система уравнений (1) в матричном виде для системы взаимоувязанных тел. В данном случае в матрице перемещений учтены накладываемые на тела связи и граничные условия, а также заданы известные внешние усилия. Полученная основная система уравнений решается шагами по времени с применением метода прямого интегрирования.

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

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

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

Методы решения задач на собственные значения и методы решения уравнения (1) шагами по времени приведены в работе [2]. Более подробное описание метода формирования матриц уравнения (1) на основании усеченного ряда собственных частот податливого тела можно найти в [1].

2. ОСОБЕННОСТИ ПОСТРОЕНИЯ МОДЕЛЕЙ ПРИ ПРОЕКТИРОВОЧНОМ РАСЧЕТЕ

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

На рис. 2 приведена одна из аэромеханических схем вертолета и его математическая модель, которая построена в расчетном комплексе ADAMS. Схема вертолета рассматривается

Рис. 2. Трехмерная математическая модель аэромеханической схемы вертолетов

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

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

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

о напряженно-деформированном состоянии элементов механической системы из неоднородных материалов;

о кинематике механической системы и величинах усилий в шарнирных соединениях и опорах;

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

о законах движения входных и выходных звеньев в стационарных и нестационарных режимах работы.

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

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

Алгоритмы решения нестационарных задач аэродинамического расчета изложены в работах Белоцерковского С. М. За прошедшие годы в ЦАГИ проведены работы по определению аэродинамической нагрузки несущего винта (например, [3], [4]). Расчетные методы направлены на описа-

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

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

Данная модель позволяет провести кинематический и силовой анализ системы управления несущим винтом для заданного режима полета. На рис. 4 изображено изменение нагрузки по времени на бустере циклического шага. Отрыв вертолета от земли и переход в набор высоты после 22 с. В данном случае расчет предназначен для моделирования типового полета по прямой.

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

Время, с

Рис. 3. Изменение подъемной силы по времени вдоль лопасти несущего винта

Рис. 4. Изменение нагрузки по времени на бустере циклического шага об/мин

Время, с

Рис. 5. Изменение количества оборотов при разгоне несущего винта

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

Кинематический и силовой анализ отдельных агрегатов и звеньев систем, а также некоторые особенности определения их напряженно-деформированного состояния рассмотрены в работе [5].

3. О МОДЕЛИРОВАНИИ СТЕНДОВЫХ ИСПЫТАНИЙ

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

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

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

Нагрузка, Н

Нагрузка, II эооо

¡Г " "

Т/ ¡1______.

■'' !/

1;----- г. ' щ;

Ч'Г~

у/

Время. I

Время, с

Рис. 6. Имитация стендовых испытаний:

а — трехмерная модель испытательного стенда несущего винта; б, в — изменение подъемной силы по времени вдоль лопасти при отличном от нулевого положения циклического или общего шага соответственно; г — движение законцовки лопасти при заданном положении циклического шага; д — изменение нагрузки по времени на бустере циклического шага (автомат перекоса изменил положение от начального по общему и циклическому шагу)

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

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

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

Балансировочный груз

а)

Гибкий торсион

Нагрузка, Н

Шй 1»

Нагрузка, П

I

и

я

■ »

Перемещения

1

Нагрузка | «

....— ¿лЛЛ 1! и

1;

¡а*

/Ц

/У/

// .У

-

Рис. 7. Имитация стендовых испытаний с учетом деформаций элементов конструкции:

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

в момент флаттера в поле центробежной нагрузки

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

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

Расчет (см. рис. 7) является дополнением к имитации стендовых испытаний, приведенной на рис. 6. В данном случае возможно получить перемещения (рис. 7, б) и изменение аэродинамической силы в одной из контрольных точек модели при набегающем на лобовую кромку прямолинейном равномерно нарастающем потоке. Возможно исследование изменений флаттерных характеристик под влиянием центробежных нагрузок. К примеру, на рис. 7, г показаны перемещения концевой части лопасти (в точке 9, см. рис. 3) по вертикали при расчете на флаттер с учетом влияния центробежной нагрузки. Так же, как и предшествующий расчет (см. рис. 6), данная математическая модель позволяет оценить деформации гибкого торсиона втулки несущего винта в плоскости вращения винта при изменении скорости вращения вала несущего винта. В соответствии с оценкой изменения подъемной силы по времени вдоль лопасти при расчете флаттера с учетом влияния центробежной нагрузки (рис. 7, в) возможна оценка законов изменения деформаций торсиона и уровня его напряженного состояния по заданной теории прочности в сочетании с изменением нагрузки в шарнирах системы управления или узле крепления гибкого торсиона (например, к валу несущего винта, рис. 7, д). На начальной стадии расчета перемещения (рис. 7, г) вызваны отклонением лопасти от горизонтального положения под действием сил гравитации, колебания продолжаются до нарастания подъемной силы по достижении лопастью заданных оборотов. При увеличении угловой скорости на лопасти проявляются предфлаттерные биения (на рис. 7, г это отрезок времени от 6 до 10 с), которые при дальнейшем увеличении оборотов переходят в нарастающие колебания (см. рис. 7, г, расходящиеся колебания после 10 с). Изменения нагрузок, соответствующие колебаниям лопасти, показаны на рис. 7, д.

Модель имеет иные варианты исполнения, например, с учетом деформаций всей лопасти в целом. Иными словами, диапазон возможностей расчетчика для исследований достаточно высок.

4. О ЧИСЛЕННОМ МОДЕЛИРОВАНИИ АВТОМАТИЧЕСКОГО УПРАВЛЕНИЯ

МЕХАНИЧЕСКОЙ СИСТЕМЫ

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

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

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

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

1) в заданной системе автоматического управления найти и оценить переходные процессы — это задача анализа системы автоматического управления;

2) по заданным переходным процессам и основным показателям разработать систему автоматического управления — это задача синтеза системы автоматического управления.

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

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

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

ВЫВОДЫ

Разработан комплекс математических моделей, который дополняет проектировочные расчеты и дает возможность моделирования некоторых режимов стендовых и летных испытаний, а также провести уточненный кинематический и силовой анализ. Математические модели разработаны на базе расчетного комплекса ADAMS с дополнительным программированием функций определения аэродинамических нагрузок и законов управления звеньями механической системы. Реализация полученных моделей возможна и на расчетных комплексах NASTRAN и ANSYS при дополнительном их программировании. Геометрические построения и расчет массовых характеристик звеньев механической системы проводились с применением комплекса твердотельного моделирования SolidWorks. Решение задач на собственные значения для учета податливости некоторых звеньев механической системы проводилось с применением NASTRAN и ANSYS с последующим использованием результатов расчета в динамической модели ADAMS.

Разработанный комплекс математических моделей позволяет наиболее полно использовать возможности численного эксперимента на ранних стадиях проектирования механических систем летательных аппаратов. Приведенные в работе модели являются дополнением к расчетным моделям, разработанным без применения NASTRAN и ANSYS [5], и позволяют увеличить мощность проектировочных расчетов в сочетании с возможностями расчетных комплексов зарубежной и отечественной разработки.

ЛИТЕРАТУРА

1. Craig R. R. and Bampton M. C. C. Coupling of substructures for dynamics analyses // AIAA J. 1968, 6(7):1313 — 1319.

2. Бате К., Вилсон Е.Численные методы анализа и метод конечных элементов. Пер. с англ. А. С. Алексеева и др. Под ред. А. Ф. Смирнова. — М.: Стройиздат, 1982.

3. Головкин М. А., Тарасов Н. Н. Аэродинамические характеристики моделей воздушных винтов вертикально взлетающего самолета в широком диапазоне углов атаки // Ученые записки ЦАГИ. 2009. Т. XL, № 4, с. 29 — 40.

4. ЛеонтьевВ. А. Метод решения уравнений движения упругих лопастей вертолетных винтов в общем случае движения // Ученые записки ЦАГИ. 2010. Т. XLI, № 5, с. 67 — 79.

5. Голованов А. И., Касумов Е. В., Шувалов В. А. О методике численных экспериментов в проектировочных расчетах механических систем вертолета // Ученые записки ЦАГИ. 2010. Т. XLI, № 4, с. 86 — 104.

Рукопись поступила 27/I2012 г.