CC BY
106
Том ХЫ
УЧЕНЫЕ ЗАПИСКИ ЦАГИ
2010
№ 4
УДК 629.735.45.064
о методике численных экспериментов в проектировочных расчетах механических систем вертолета
А. И. ГОЛОВАНОВ, Е. В. КАСУМОВ, В. А. ШУВАЛОВ
Рассматривается методика проектировочного расчета механических систем. Предлагаемая методика построена на объединении нескольких групп задач численного моделирования при решении задач статики и движения механических систем под воздействием внешних сил. Приводятся примеры расчетов конструкций из композиционных материалов, механизмов системы управления легкого вертолета. Рассматриваются возможности применения решения нелинейных задач на ранних стадиях проектирования летательных аппаратов.
Построение алгоритмов решения задач проектирования основано на смешанном использовании программных расчетных комплексов Российской и зарубежной разработки.
Ключевые слова: численный эксперимент, математическая модель, проектирование, вертолет, композит, композиционный материал, механическая система, кинематика, силовой анализ, метод конечных элементов, деформация, напряжения.
Проектирование и испытания современных механических систем носят комплексный характер. Испытания на натурных объектах широко совмещаются с методами математического моделирования с использованием ЭВМ. Это привело к возможности развития методов планирования натурных испытаний, которые обеспечивают большую точность получаемой информации.
Целью настоящей работы является поэтапное создание комплекса нелинейных математических моделей проектировочного и поверочного расчетов механических систем вертолетов. Математические модели разрабатываются на ранних стадиях проектирования для моделирования некоторых режимов до построения стендов и опытных образцов летательного аппарата. По мере накопления материалов стендовых и летных испытаний предполагается возможность уточнения математических моделей для последующих поверочных расчетов.
Говоря о проектировании механических систем в целом, необходимо сказать, что в настоящее время развитие получают машины автоматического или полуавтоматического управления, обладающие реакцией на неощутимую человеком информацию, например на ультразвук, вибрации, электромагнитные и тепловые поля и т. п. Отличительной чертой автоматических механических систем является наличие логических элементов автоматической подстройки (электронного, пневматического, гидравлического и механического типов) управляемых объектов при постоянно изменяющихся параметрах внешних условий для качественного выполнения требуемых режимов работы.
При проектировании машин особое внимание уделяется задачам определения вибрационных характеристик. С одной стороны, это вопросы борьбы с вибрациями путем создания виброустойчивых механизмов, с другой стороны, это исследование резонансных характеристик вибраций для выполнения заданных функций механизма, обладающего требуемыми кинематическими характеристиками. Вибрации искажают основное движение элементов машин, механизмов и систем управления по предписанным кинематическим законам, порождают неустойчивость заданного закона движения и часто приводят к отказу всей системы.
В настоящей работе построение алгоритмов решения задач проектирования основано на совместном использовании программных расчетных комплексов российской и зарубежной разра-
ботки (NASTRAN, ANSYS, ADAMS) с целью оценки их положительных и отрицательных сторон, а также возможности их эффективного сочетания при реализации решения нелинейных задач механики, аэромеханики и аэродинамики.
Иными словами, первоначально подбирался набор численных методов решения необходимых систем уравнений (набор конечных элементов [6, 7], алгоритмы решения задач на собственные значения [1], метод определения аэродинамических нагрузок [2, 3, 9] и т. д.). Следующим этапом являлась проработка последовательного применения алгоритмов определения нагрузок (аэродинамических, температурных и т. д.), расчета динамики системы тел в единой математической модели. Проведен анализ точности решений единого взаимно увязанного комплекса задач вычислительной математики по определению нагрузок различного вида, кинематики и силового анализа, методов расчета деформаций тел в составе механической системы. В заключение, на основании подобранного по точности решения математического аппарата сформулированы алгоритмы проектировочного расчета реальных механических систем различного типа.
При решении нелинейных задач вычислительной математики была рассмотрена система управления лопасти несущего винта легкого вертолета из композиционных материалов.
Принцип действия несущей системы вертолета и ее взаимодействие с окружающей средой определяют широкий спектр задач математического моделирования. С одной стороны, вертикальное и поступательное движения вертолета обеспечиваются за счет изменения угла атаки лопасти несущего винта. С другой стороны, колебания несущей системы должны сочетаться с частотными характеристиками конструкции вертолета в целом и не вызывать резонансных явлений на различных режимах полета, при взлете и посадке.
Применение композиционных материалов в конструкции лопастей и звеньев механизмов управления приводит к необходимости учета изменения механических свойств материалов под воздействием температуры и влажности в сочетании с изменением давления по высоте полета в поле переменных центробежных и аэродинамических нагрузок. Программирование входных звеньев механизмов системы управления в данном случае определяет запасы управления летательным аппаратом, его маневренность и устойчивость.
Особенности формирования исходных данных. Система управления лопасти несущего винта с бесшарнирной втулкой вертолета (рис. 1, б) по своим конструктивным особенностям отражает основные направления развития современных механических систем и позволяет наиболее полно апробировать алгоритмы расчетов [5]. Она характеризуется тем, что:
применяются композиционные материалы;
кинематика системы управления основана на деформациях гибких звеньев (торсионов) втулки несущего винта;
входные звенья системы управления (гидравлические бустеры) управляются электронными блоками, которые работают по заданным законам управления в зависимости от температуры и давления окружающей среды и углов отклонения рычагов управления.
Лопасти и втулка несущего винта вертолета изготовлены из стеклопластика (рис. 1, г). Втулка состоит из (рис. 1, а): 1 — корпуса втулки и деталей, соединяющих торсионы втулки с валом винта; 2 — четырех кожухов управления для управления шагом лопастей; 3 — четырех металлических узлов крепления лопасти в виде переходников с проушинами для торсиона и лопасти; 4 — двух упругих композиционных балок, в каждом рукаве которых имеется гибкое звено — торсион.
Заготовку упругой балки торсиона (рис. 1) собирают из предварительно спрессованных пакетов стеклоткани и листов резины, чередующихся по высоте. Затем методом горячего прессования изготавливают упругую балку с припуском на каждую сторону. Окончательные размеры получают механической обработкой. Угол укладки слоев стеклопластика к продольной оси торсиона ноль и ±45°.
Автомат перекоса (рис. 1, в) выполнен по трехточечной схеме и имеет центральный сферический подшипник с оргалоновым покрытием. Углы наклона кольца управления задаются сочетанием длин штоков бустеров для каждого режима полета.
Переход из одного режима полета вертолета в другой задается углом наклона кольца управления автомата перекоса вследствие изменения длины штоков гидравлических бустеров. Бустеры управляются сигналами электронной системы управления на основе решения линеризованной
Рис. 1. Конструкция упругой балки втулки и автомата перекоса системы управления несущим винтом: 1 — корпус втулки; 2 — кожухи управления; 3 — металлические узлы крепления лопастей; 4 — рукав упругой балки (торсион)
системы уравнений, которая отражает законы управления аэромеханической схемой вертолета. Входными параметрами для электронной системы являются углы отклонения рычагов управления и переменные параметры внешних условий (давление и температура).
В идеале, разрабатываемый комплекс математических моделей должен давать максимально полное представление:
о напряженно-деформированном состоянии элементов механической системы из неоднородных материалов;
о кинематике механической системы и величинах усилий в шарнирных соединениях и опорах;
о частотных характеристиках системы в целом при изменении внешних параметров и нагрузок;
о законах движения входных и выходных звеньев в стационарных и нестационарных режимах работы.
Математическое моделирование автомата перекоса проведено в предположении абсолютной жесткости его звеньев. Основное дифференциальное уравнение движения выражается в матричной форме для системы взаимоувязанных тел с учетом условий их сопряжения, граничных условий и системы внешних нагрузок. Расчет позволяет рассматривать возможности влияния сил трения и наличия люфтов в шарнирных соединениях на нагружение звеньев при различных режимах полета вертолета. Для примера на рис. 2, а приведена изменяемая во времени расчетная нагрузка на автомат перекоса несущего винта вертолета от действия аэродинамических сил на заданном режиме полета. При необходимости математическая модель позволяет применять данные о внешних нагрузках из материалов натурных испытаний. В данном случае аэродинамическая нагрузка определена заранее для заданного режима полета по методу присоединенных вихрей без учета влияния поверхности конструкции вертолета на лопасть несущего винта.
В качестве примера на рис. 2, б продемонстрирован расчет изменения во времени значения реакции в опоре одного из бустеров управления автоматом перекоса в проекциях на оси координат. Данный расчет проведен без учета скорости перемещения рычагов управления. По этой причине углы установки звеньев автомата перекоса для заданного режима полета устанавливаются в период 0.1 с (это видно по характерному скачку на графике изменения расчетной величины усилия), затем определяется искомая величина усилия. При отработке решения с учетом разных скоростей движения рычагов управления кривая усилия будет сглаживаться в соответствии со скоростью движения штоков бустеров управления.
Рис. 2. Математическая модель системы управления несущим винтом без учета податливости гибких звеньев
Последующие расчеты на статическое нагружение для характерных режимов позволяют оценить напряженно-деформированное состояние звеньев механизма.
Рассмотренная динамическая модель не позволяет оценить кинематику системы управления несущим винтом с учетом деформаций втулки несущего винта и лопастей. В этом случае моделирование геометрии гибкого звена проводится методом конечных элементов для предварительного определения его частотных характеристик при возможных вариантах граничных условий. Последующее решение задачи движения взаимоувязанных звеньев на основе их частотных характеристик дает более полное представление о кинематике механизма в целом и его динамическом отклике на действие переменной нагрузки.
Рассмотренный вариант системы управления несущим винтом не является единственным [4, 8]. Возможно моделирование иных кинематических схем.
Математическое моделирование механической системы. Для решения подобной задачи с учетом податливости звеньев (рис. 3) необходимо применять специальные математические подходы. Один из них изложен в работе [11], в которой уравнения движения для гибкого тела получены из уравнения Лагранжа в форме:
(дЬ ^
Ж
дё
дь с¥
я-е=о,
(1)
¥ = о,
(2)
где Ь — Лагранжиан; ¥ — функция разложения энергии; ^ — граничные условия; X — множители Лагранжа; ё — обобщенные координаты; е — обобщенные приложенные усилия.
Ь = Т - V,
(3)
где Т, V — кинетическая и потенциальная энергии соответственно.
Для кинетической энергии в обобщенной матрице масс и обобщенных координатах системы справедливо выражение:
1 ;
Т = - ёТ м (ё)ё.
(4)
Рис. 3. Вектор в глобальной системе О, который описывает положение точки Р (изменение положения точки Р в точку Р' на искаженном теле В определяется относительно местной системы координат)
Матрица масс М (Д) представлена матричными блоками 3 х 3:
М ф =
Мп Мг Мт
МТ МГ1
Мг
мТ мТ м
Ш гт тт
(5)
где нижние индексы I, г и т обозначают перемещения, повороты и модальные значения соответственно. Модальные значения определяются методом модальной суперпозиции, они рассмотрены ниже. Матрица масс М выражается в девяти инвариантах инерции:
М„ = ,
м г =-А
А2 + JЪ]
В,
М т = А33
М гг = В1
37 -
А 8+А 8Т
§,- А
3 "3^3
(6)
В,
М гт = В1
А4 + Чз
М тт = ^6.
Здесь А и В — матрицы преобразования системы координат, ч — модальные координаты. Очевидна явная зависимость матрицы масс от модальных координат.
Инварианты инерции вычислены для N узлов конечно-элементной модели, основанной на информации о массе каждого узла тр.
Потенциальная энергия выражается в обобщенных координатах:
у = Уё 00+2 Щ. (7)
В выражении потенциальной энергии во втором слагаемом выражаем эластичность тела, К — обобщенная матрица жесткости, которая постоянна. Только модальные координаты ч вносят вклад в упругую энергию (см. (14)).
Поэтому форма матрицы К имеет вид:
к оо=
Кп КГг Кт "0 0 0 "
Кгг К гт 0 0 0
Кт Кт гт К тт _ 0 0 К тт _
где Ктт является обобщенной матрицей жесткости структурного компонента относительно модальных координат q.
В первом слагаемом (7) величина V — это гравитационная часть потенциальной энергии:
Vg =\ргр • gdV = {р[х + А( + Ф(Р)q)gdV,
(9)
где g — вектор ускорения свободного падения.
Результирующая гравитационная сила имеет вид:
6У£ дё
дА
дё
¡рсИУ
V
{р( +Ф(P)q) (¿V |рФт(P)dV
(10)
Силы демпфирования зависят от обобщенных модальных скоростей и выражены в форме:
(11)
F = 2 Щ,
где матрица Б — симметрична, постоянна и содержит коэффициенты демпфирования .
В случае ортогональных модальных форм матрица демпфирования может быть эффективно определена через диагональную матрицу модальных отношений демпфирования е1 (см. [11]).
Обобщенная сила О состоит из обобщенной приложенной силы, обобщенного момента и обобщенной модальной силы в виде:
О
^т QR Яы
(12)
Обобщенная модальная сила — модальная сила на теле, приложенная в точке или указывающая момент, полученный в проекции силы на модальную форму.
Для кинематического и силового анализа механической системы с учетом податливости ее звеньев из композиционных материалов дифференциальное уравнение движения выражается в обобщенных координатах в следующей форме:
+ М ё - 2
дМ
дё ё
ё + Кё + +Бё-
дё
(13)
где ё, ё, ё — обобщенные координаты податливого тела и их производные по времени; М — матрица масс податливого тела (5); М — производная по времени матрицы масс податливого
тела; -^М — частная производная матрицы масс гибкого тела относительно обобщенных координат, где М — число мод; К — обобщенная матрица жесткости; — обобщенная гравитационная сила; Б — модальная матрица демпфирования.
Основным предположением для описания поведения гибкого тела является то, что рассматриваются малые линейные деформации тела относительно местной системы координат, в то время как местная система координат, в свою очередь, испытывает большое нелинейное перемещение в глобальной системе координат.
Линейные деформации тела в локальной системе координат рассматриваются как линейные деформации узлов элементов конструкции, приближенно представленные линейной комбинацией конечного числа векторов формы ф:
м
= 1ФЧ. (14)
1=1
Здесь М — номер формы, ч — фактор масштаба (амплитуды) модальной координаты. Сложная форма строится как линейная комбинация простых форм, что демонстрируется на рисунке ниже:
Данный метод носит название метода модальной суперпозиции. Его основное предположение — компоненты деформаций с большим рядом модальных значений могут быть рассмотрены при сокращенном числе мод. Сокращенное число мод называют модальным усечением или модальным основанием. Иными словами, для решения задачи моделирования механической системы спектр собственных векторов рассматривается в усеченном виде при различных вариантах граничных условий (набор граничных условий).
В методе модальной суперпозиции применяется понятие уровень динамического содержания гибкого тела — это набор частотных характеристик гибкого тела с учетом особенностей изменения граничных условий в механической системе. Чем полнее определено модальное основание, тем точнее гибкое тело при решении основного дифференциального уравнения движения отражает частотные характеристики и деформации реального звена механической системы.
Уравнение (14) в матричной форме имеет вид:
и = ФЯ, (15)
где ц — вектор модальных координат и моды в виде модальной матрицы-столбца Ф. После модального усечения Ф становится прямоугольной матрицей.
Модальная матрица Ф является преобразованием от маленького набора модальных координат ч к большему набору физических координат и.
Модальное усечение разделено на две части — граничную и основную (рис. 4).
Граничная часть модального усечения — это сочетания граничных условий, которые полностью охватывают возможные перемещения границ податливого тела, чс = ив.
Основная часть модального усечения — это модальные формы для каждого рассматриваемого варианта граничных условий. Расчетчик имеет возможность выбора сколь угодно большего количества модальных форм с целью повышения точности расчета.
а
Ь
Рис. 4:
а — набор граничных условий, состоящий из двух типов ограничения левого конца луча (луч имеет две точки приложения на обоих концах, левая форма соответствует перемещениям одного из узлов, правая — вращению в узлах на концах луча); б — две основные модальные формы для луча при граничных условиях с закреплением его на обоих концах
Отношение между физическими модальными формами и формами по методу суперпозиции выражается следующим образом:
I
Фс
0
Фт
(16)
где пв — граничные модальные усечения; и1 — основные модальные усечения; I, 0 — идентичность и нулевая матрица соответственно; Ф^ — физические смещения основных модальных усечений в наборе граничных условий; — физические смещения основных модальных усечений в нормальных формах; дс — граничная часть модального усечения; ^ — основная часть модального усечения.
Обобщенные матрицы жесткости и масс выражаются в данном методе через полученное модальное основание.
Матрица жесткости имеет вид:
К = Ф7 КФ
I 0 " Т КВВ КВ1 " I 0 "
Ф1с Фт _ _ К1В Кц _ Ф1с Фж _
К
сс
0
0 К
NN
(17)
Матрица масс:
М = Ф7 МФ
I 0 " Т 'МВВ МВ1" " I 0 "
Ф1с Фж _ _ М1В Ми _ _Ф1с Фж _
М М,
сс
CN
М М
Ж
NN.
(18)
где нижние индексы I, В, N и с обозначают основные, граничные модальные усечения и набор граничных условий. Матрицы М, К — это обобщенная матрица масс и матрица жесткости соответственно.
Уравнения (17), (18) имеют несколько особенностей: Мж, Кж — диагональные матрицы; К — блочная диагональная матрица, которая не имеет никакой связи со способами ограничения и границами нормальных мод в модальных усечениях; М — диагональная неблочная матрица, потому что имеется связь в части инерции между граничными способами и нормальными модами в модальных усечениях.
Метод модальной суперпозиции — это метод для создания модального основания, достаточного для описания гибкости тела при решении основного уравнения движения механической системы. Однако изначальные данные для метода имеют ряд недостатков:
1. Необходимо учитывать, что в системе имеют место шесть перемещений как твердого тела, которые должны быть устранены перед решением основного дифференциального уравнения движения.
2. Набор граничных условий в методе суперпозиции — результат статической конденсации. Следовательно, эти моды не отражают динамического поведения в полной мере для гибкого тела.
3. Набор граничных условий невозможно рассмотреть иначе, в противном случае мы получили бы ограничения на систему уравнений в целом.
и
Эти проблемы в методе решаются с применением специальных математических действий.
Метод модальной суперпозиции — это набор граничных условий и частотных характеристик к ним. По этой причине перед решением основного дифференциального уравнения движения непосредственно для гибкого тела решается задача:
Кя = Шц. (19)
Решив ее, мы получаем собственные векторы, которые приводим с помощью матрицы преобразования N к эквивалентному, ортогональному основанию с модальными координатами д*, для метода суперпозиций:
N4* = я. (20)
Тогда метод суперпозиций выражается формулой:
мм м
и = !Ф я =Х>» N4* =£*У, (21)
г=1 г=1 г=1
где ф* — это собственные векторы исходной системы, они являются собственными векторами
метода суперпозиций для связи с естественными частотами системы.
Описать физический смысл этих мод трудно, но можно сказать следующее. Установленные границы физически реальных мод заменены приближенно собственными векторами податливого тела. Это приближение, так как оно основано только на модальном усечении. Из этих мод шесть — моды абсолютно жесткого тела. Граничные условия заменяются границами собственных векторов.
Более подробное описание применения метода суперпозиций можно рассмотреть в работе [11].
Особенности применения метода конечных элементов в моделировании кинематики механических систем [10]. Как отмечалось выше, К — блочная диагональная матрица жесткости (17), которая не имеет никакой связи со способами ограничения и границами нормальных мод в модальных усечениях. Это обстоятельство позволяет разрабатывать различные варианты формулировок матриц жесткости конечных элементов (45) с целью повышения точности решения в отрыве от метода решения основного дифференциального уравнения движения.
В геометрически нелинейной постановке прогибы элементов конструкции соизмеримы с толщиной или превышают ее в несколько раз. С физической точки зрения речь идет о моделировании процесса состояний, каждое из которых реализуется при определенном значении внешней нагрузки. В общем случае на каждом шаге нагружения может меняться и структура распределения внешних сил.
С учетом особенности решаемой задачи и перспектив ее дальнейшего развития необходима разработка специального трехмерного конечного элемента, имеющего слоистую структуру по толщине с различными механическими характеристиками каждого слоя.
В этом случае построение матрицы жесткости имеет свои особенности, так как общий интеграл потенциальной энергии деформации представляется в виде суммы соответствующих значений по слоям. Достоинство подобной модели состоит в более точном представлении напряженно-деформированного состояния элемента (по слоям), и это позволяет избавиться от погрешности, связанной с осреднением физико-механических свойств пакета слоев. Точность решения, получаемого в результате использования такого подхода, зависит от степени аппроксимации поля перемещения в направлении поперек слоев. Если свойства слоев сильно отличаются, то использование простых аппроксимаций в пределах элемента может вносить значительную погрешность. В этом случае необходимо либо применять полиномы высокой степени в каждом элементе и использовать один элемент по толщине, либо задавать несколько элементов по толщине, построенных на основе простых аппроксимаций. С точки зрения общего числа степеней свободы оба подхода отличаются немного, хотя в точности распределения перемещений, деформаций и напряжений отличие наблюдается.
Рис. 5
Структурно предлагаемый конечный элемент (рис. 5, а) представляет собой искривленный параллелепипед с линейчатой поверхностью в поперечных срезах, состоящий из набора слоев по толщине, каждый из которых является ортотропным материалом с осями ортотропии а, в, у.
Пошаговая схема решения геометрически нелинейных задач теории упругости (рис. 5, в) подразумевает наличие некоторой глобальной декартовой системы координат хь Х2, Х3 с ортами е1, е2, е3, в которой рассматривается процесс деформирования тела, нагружаемого заданными массовыми Q и поверхностными Я силами. Это нагружение разделяется на некоторое количество шагов и для каждого из них определяется уровень действующих сил: 0(т)
п(т) и Я ' —
(т)
радиус-вектор материальной точки на ка-
значение внешних сил на т-шаге нагружения; г
-40)
ждом шаге; п ' — положение рассматриваемой точки в начальном состоянии.
Для описания деформирования всего тела будем использовать Лагранжев подход. Для этого введена «вмороженная» система координат ^з , относительно которой положение каждой
точки при движении неизменно, так как «деформируется» сама система этих координат. Уравнение равенства работы имеет вид:
ШХ 4т Цт) * а _ Ж (т О +|Ц Я(т)ЬШ8,
(22)
Оп ', з
где — часть поверхности, где заданы поверхностные силы.
Для получения разрешающего уравнения вводится в рассмотрение приращение деформаций:
ДЕ(т) _ Е(т) - Е(т-1) _ 1 ~ Еи Еи ~ 2
дДЖ(т) дДЖ,(т)
дДГк(т) здт
(т-1)
к
ддж(т) ддж:
'к
к
дх,
.(0)
дх
.(0)
дх
.(о)
дх.
(о)
дх,
(о)
(т-1)
(т) дД#(т)
дД^Г> дДЖ
дх
'к .(0)
дх'
к (0)
(23)
к
В этих обозначениях вариационное уравнение имеет вид:
ЫЕ( 4т-1)+л^л«^=Ы й (т)5 ^+ц я(т ЫёБ.
Оо >, 1
В отличие от уравнения (22), уравнение (24) допускает линеаризацию в предположении, что градиент приращений существенно меньше единицы, т. е.
дЛУУ(т)
4°)
«1.
(25)
Для этого определим линейную часть от ЛЕ1 ' в виде:
( ) 1 17 дШ (т-1)^ЛЛ
А ~(т ) 1 х-1
Л4т) = 21
5
д,к я, дх.
(о)
(т) дЛУ(т
дх
(о)
1к
дх(,0)
1 /
дх,
.(о)
(26)
С учетом этого обозначения, после пренебрежения малыми второго порядка получим линейное уравнение:
(
ЖЕ
оп 1,1
Л4т )5в(т) +1 )5
1 ^ дЛУкт) дЛ1гкт)
2 Е-
дХ((0) дх(,0)
Л
й 0 =
(27)
: й(т ]5Ш О + Ц Я (т)5ШБ - Ц|Е 1 )5в ((т )й О.
Оо 1,1
Второй вариант предполагает в качестве неизвестных функций вектор V(т), который является вектором приращения перемещений при переходе из (т - 1)-го состояния в т-е, т. е.
V(т )=Л\¥(т). (28)
Вариационное уравнение принципа возможных перемещений, аналогичное (22), (24), имеет вид:
Ш Е(4т-1)+Л4т))5е(т)й0= й(т)5^0+ Ц Я(т)5ЙБ. (29)
0(т-1)] 0(т-1) Б™-1)
В предположении справедливости (26) с учетом (28) уравнение (32) тоже допускает линеаризацию. Если ввести «классические» линейную и нелинейную деформации:
>) = 1 ^ = 2
(
д^т) д^т)
\
(т-1) дх(т-1)
V 1 1
п(т)
% "2Г дх(т-1)дх(т-1), то после линеаризации (29) полученное линейное уравнение имеет вид:
Ш Е(ЛБ<т)5в(т)+а^)5п(т))й0 =
о(:-1) /, 1
= ^ й(т)5^0+ Ц Я(т)5ШБ- Ест(/т-1)5в!7т)й0.
(3о)
(31)
(т-1)
(т-1)
(т-1) 1, ]
Отличие уравнения (22) от аналогичного (27) состоит в более простом соотношении для линейных деформаций. Однако это упрощение сопровождается дополнительным вычислением напряжений с/т) для перехода к следующему шагу нагружения. Для этого необходимо воспользоваться соотношениями, связывающими тензоры напряжений Коши и Кирхгофа, которые в нашем случае имеют вид:
/ т) _.
1
с;,- 7 _
- /т) дх/т) /
3«/ дХ) /с/т-1)
./ т-1) к!
к ,1 дх^" 1) дх ,
-ДЯ
/т)
к!
),
(33)
где
3
дх1т) дх/т) дх3т)
дх1/т-1) дх2т-1) дх}т-1)
дх(т) дх/т) дх3т)
дх/т-1) дх2т-1) дх/т-1)
дх1т) дх/т) дх3т)
дх3т-1) дх3т-1) дх3т-1)
(34)
В соответствии с общей схемой метода конечных элементов исходное тело представляется в виде совокупности отдельных конечных элементов, в каждом из которых для неизвестных функций определяются самостоятельные аппроксимации. Так как в процессе пошагового нагру-жения геометрия элемента претерпевает изменения, то для определения радиуса-вектора геометрии используется тоже аппроксимация. Очевидно, что наиболее удобно воспользоваться аппроксимацией одного порядка как для геометрии, так и для перемещений, т.е. использовать конечные элементы изопараметрического типа.
Используемые аппроксимации имеют вид:
{ У )}_1
{ / у )}_1
«(т-1)
1, п
«(т-1)
2, п
х/т-1) с3, п
у( ™ )
1, п
у(т)
2, п
уН
3, п
*п / У )_!{хПт-1)}^п / у )
(35)
N / ] )_1{пт }Ып / ] )
где Ып /у ) — функции формы относительно безразмерных координат Еу /-1 < Еу < +1), которые
можно считать лагранжевыми координатами конкретного элемента, т. е. они совпадают с теми лагранжевыми координатами Еу, которые были определены ранее. Связь между произ-
водными по Е через производные по х1
■(т)
определяется с помощью матрицы Якоби
^ ]_1
дЕх
дЫп дЕ2
дЫп
дЕ3
{х/т-1)}7
(36)
п
которая определяет связь:
д^
Обратная матрица Якоби имеет вид:
д
дх5т-1)
дх(т-1)
-1
д?
Таким образом, производные:
дNn
дХ((т-1)
-1
N.
д?,
Введем векторные аналоги компонент тензоров деформаций и напряжений: ~(т)}Т = |р(к) р(к) р(к) р(к)+р(к) р(к)+р(к) р(к)+р(к)}
, р22 , р33 , р12 +р21, р23 +р32 , р31 + р13
Р"}' ={р(1
{(т-1)} ={_(т-1) .(т-1) р(т-1) _(т-1) .(т-1) _(т-1) Г ) 11 22 33 ,12 , 23 , ст31
{лб (т)}т = {лб1(1:), лб«, л4:), л4:), лб«, лб(3
Аппроксимации для линейных компонент деформаций имеют вид:
К »ЬЕЫдК >},
Ч )(?1)
где для матриц
получаются выражения:
В
(п)
дNn 0 0 дNn 0 д^
дх(т-1) дх2т-1) дх3т-1)
Т 0 дNn 0 дNn дNn 0
дх(т-1) дх(т-1) дх3т-1)
0 0 дNn 0 д^ д^
дх3т-1) дх(т-1) дх}т-1)
(38)
(39)
(40)
(41)
(42)
Для дальнейших построений необходимо определить связи между деформациями и напряжениями (соотношения упругости). В рассматриваемом варианте шагового решения эти соотношения должны связывать деформации р3 с приращениями напряжений Л8(т). В первом приближении их можно считать линейными, т. е. деформирование на шаге нагружения является геометрически и физически линейным. Физическую нелинейность в этой постановке можно учесть самостоятельным определением соотношений упругости на шаге нагружения в каждой квадра-
турной точке с учетом накопленных деформаций и напряжений. Формально эту зависимость можно записать следующим образом:
{as( " )(i, )}=[d( " у )]{s< - >(5,)}.
(43)
Конкретный вид матриц [.О] зависит от свойств материала.
Подставляя выражения (42), (43) в (32), можно записать матричный аналог этого уравнения. В частности, первый интеграл на уровне одного конечного элемента приводит к матрице жесткости и может быть записан как:
{as(-i} {ssl-)})
(--1) i, j
K
(—)
V,
(—)
(44)
где
K
— блоки матрицы жесткости рассматриваемого элемента на —-м шаге нагружения.
Для их вычисления используется численное интегрирование по некоторой квадратурной формуле, вид и порядок которой зависит от формы конечного элемента и порядка аппроксимации в функциях формы. Структурно для них можно записать следующее выражение:
K
(—)
Jf J[Bi(5 j )f [ D(—)(5 j )][B„ (5 j )] det [J (5 j )] 5
Xb (5<-))]Г ГЫ'"^' ))][b „ (5<-))] det [J (5())
Ю
t>
(45)
где 5 j — координаты квадратурных точек, rat — квадратурные множители.
Для определения второго слагаемого в левой части уравнения (32) составим выражение вариаций нелинейных компонент деформации (31) с учетом аппроксимаций (35). В частности:
где
(46)
i,n V 5V(—{vn—)}',
dNt dNn
cX<— -1) ck}—-1)
dNt dNn
dx2— -1) dx2—-1)
dNt dNn
dx3— -1) dx3—-1)
dNt N dNt dNn
cxf--1) 1)1 -1) dxi--1)
dNt dNn dNt dNn
Йг2—-1) ск-- 1)1 dx3—-1) dx2—-1)
dNt dNn dNt dNn
dx3—-1) ск— 1) ' dxi--1) dx3—-1)
(47)
Q
Если записать второй интеграл в левой части (32) на уровне одного конечного элемента в
виде суммы, аналогично (44), то появятся блоки матрицы геометрической жесткости для вычисления которых используют структурную схему их вычисления типа (45):
К
у(т)
К
*(т)1 =
Г/1 {{->(? 1 )}Т { 1} г / 1)] й ?
1 0 0 0 1 0 0 0 1
(48)
:Е{4т-1)}Т { ^ ]]}-е. [ ^ $]]
и,
1 0 0 0 1 0 0 0 1
где <с
(т-1) ] —
I
значения компонент накопленных напряжений Коши в '-й квадратурной точке
рассматриваемого элемента.
Правая часть уравнения (32) может быть представлена в векторном виде для одного конечного элемента как:
Е{5^т)} {р(т]},
(49)
где введен обобщенный вектор узловых сил {р/т)} , для которого справедливо выражение:
+1
(50)
и.
{"}=///V N1 (? 1 ]{ 1 ]}-[ вД? > )]Т { > ]}} ^ Г' 1 ]] й ?1
+ // N1 ]{ ]}=
ыП-1
= ЕГN. 1']){йМ(?/>]}-[(]]]Т К'"4}]-е.Г/]]
, V -1 у1-
+ Е N (]{{(?•' ]}' ]и.
Здесь и 1 ', ? 1 ' — квадратурные множители и координаты соответствующих точек на границе с заданными граничными усилиями.
После того как определены все слагаемые на элементе, строится процедура сборки конечного элемента в единую конструкцию путем специального суммирования матриц (45), (48) и век-
торов (49) в глобальные матрицу жесткости
К
(т)
матрицу геометрической жесткости
К
(т)
и вектор узловых сил {р(т)}. Если обозначить через {у(т]} глобальный вектор узловых перемещений, то получим систему линейных алгебраических уравнений:
(
К
(т)
К
(т ]
]{у
(т Л = |Р(т ]
ГГ ]}-
которая будет дискретным аналогом вариационного уравнения (32).
Далее вводятся кинематические граничные условия (обычно за счет специальной коррек-
ции глобальных матриц
К
(т)
К
(т ]
), решается система линейных алгебраических урав-
нений, находится вектор узловых перемещений {у(т)} и определяется переход из (т - 1)-го состояния в 3-е.
Переход к следующему шагу по нагрузке требует вычисления новой геометрии и напряжений Коши. На уровне конечного элемента процедура сводится к следующим соотношениям: — вычисление новых значений координат узлов:
= {х(т-1)1
у
(т П
{х(: ]}-{А П ^{У }, — вычисление приращений напряжений Кирхгофа во всех квадратурных точках:
{Л8(т ]} = Г В(т ](?1 ]Мв п (?1 ]Шт)},
(52)
(53)
вычисление матрицы
(33) в квадратурных точках:
4=1-
дNn
дг(т-1)
Г(т) к,п '
(54)
вычисление новых компонент тензора напряжений <с по
{с(т)} по формуле:
(т) _.
с>,- _
X ^гк J1
1 (<
.(т-1) .
к1
-лб:
(т) к1
(55)
определяющих вектор <с
Таким образом, напряженное состояние описывается совокупностью напряжений в каждой квадратурной точке каждого конечного элемента и обновляется от одного шага нагружения к другому.
На рис. 5, г продемонстрировано распределение продольных напряжений с^ на внешней поверхности торсиона при статическом нагружении с заданными граничными условиями [6].
О математическом моделировании механических систем с гибкими звеньями. Применительно к расчету механизмов следует учитывать постоянное циклическое изменение по величине и направлению нагрузки на звенья. В рассматриваемом ниже примере системы управления несущим винтом вертолета подобный расчет позволяет оценить взаимодействие жестких звеньев автомата перекоса и гибкого торсиона под влиянием заданных нагрузок. При решении основного уравнения движения для механической системы нагрузки изменяются на каждом шаге интегрирования по времени в зависимости от угла атаки лопасти несущего винта.
На рис. 6 представлена математическая модель системы управления несущим винтом легкого вертолета с учетом податливости звена из композиционных материалов (торсиона) при абсолютно жестких звеньях автомата перекоса и кожуха управления для управления шагом лопасти.
На рис. 7 показана математическая модель системы управления несущим винтом легкого вертолета с учетом податливости звена из композиционных материалов (торсиона) при абсолютно жестких звеньях автомата перекоса, кожуха управления и лопасти несущего винта.
При описании матрицы жесткости податливого звена механической системы возможно применить изложенные выше соотношения конечного элемента, либо уравнения иного конечного элемента, удобного расчетчику. При решении кинематической задачи на основе математической
п
Рис. 7
модели, представленной на рис. 6 и 7, расчетчик имеет возможность оценить перемещения узлов конечно-элементной сетки податливого звена. Определение уровня напряжений для заданного момента времени требует дополнительного расчета, аналогичного изложенному выше.
Апробирование модели, приведенной на рис. 6, показало возможность расчета нагружения звеньев конструкции на различных режимах полета с учетом изменения во времени аэродинамических и центробежных сил.
Рассмотрены варианты модели с установкой несущих поверхностей и расчетом аэродинамической нагрузки по методу присоединенных вихрей (рис. 7). В данном случае кинематический и силовой анализ механической системы проводятся относительно переменной аэродинамической нагрузки, переопределяемой на каждом шаге интегрирования в зависимости от изменения углов атаки лопастей. Результаты расчетов такого варианта решения в публикации не приводятся. Данная задача требует подробного рассмотрения в отдельной работе, так как существует специфика расчета системы внешних нагрузок на несущую систему вертолета с учетом особенностей аэромеханической схемы. Требуется подробное изложение применения метода определения аэродинамической нагрузки с учетом устойчивости аэромеханической схемы летательного аппарата и введенных в математическую модель допущений относительно взаимного влияния элементов конструкции.
Разработанные модели позволяют проводить оценку некоторых режимов испытаний несущей системы вертолета на ранних стадиях проектировочных расчетов до построения натурных стендов и опытных образцов летательного аппарата.
Выводы. Проектировочные расчеты, используемые на ранних стадиях проектирования (до построения опытных образцов стендовых и летных испытаний), не позволяют достаточно полно оценить технические характеристики агрегатов проектируемой механической системы. Это приводит к дополнительным поверочным расчетам после накопления результатов стендовых и летных испытаний и дополнительным доработкам опытных образцов. Как правило, доработки приводят к увеличению веса конструкции.
В результате проведенного исследования разработан комплекс математических моделей, который дополняет проектировочные расчеты и дает возможность моделирования некоторых режимов стендовых и летных испытаний.
На ранних стадиях проектирования до построения опытных образцов комплекс нелинейных математических моделей позволяет более полно оценить:
напряженно-деформированное состояние элементов механической системы из неоднородных материалов;
кинематику механической системы и величины усилий в шарнирных соединениях и опорах;
частотные характеристики системы в целом при изменении внешних параметров и нагрузок;
законы движения входных и выходных звеньев в стационарных и нестационарных режимах работы.
Комплекс математических моделей состоит из трех групп.
Первая группа — модели кинематического и силового анализа отдельных агрегатов механической системы. Расчет проводится при условии абсолютной жесткости звеньев агрегата. Модели без учета податливости композитных звеньев (см. рис. 2) не в полной мере отражают реальную кинематику механической системы в целом и динамический отклик конструкции на действие переменной нагрузки.
Вторая группа — модели для определения напряженно-деформированного состояния и частотных характеристик отдельных звеньев агрегатов механической системы при заданных вариантах граничных условий. Расчет проводится методом конечных элементов в линейной и нелинейной постановках для заданных вариантов граничных нагрузок и случаев нагружения. Полученные модели данной группы (см. рис. 5) не в полной мере отражают напряженно-деформированное состояние и частотные характеристики отдельных звеньев, так как граничные условия определены недостаточно полно.
Третья группа — модели для проведения кинематического и силового анализа с учетом податливости звеньев из композиционных материалов (см. рис. 6, 7). Расчет проводится с применением метода модальной суперпозиции по результатам расчета частотных характеристик композитного звена методом конечных элементов для заданных вариантов граничных условий. Агрегаты механической системы в усложненных моделях рассматриваются взаимоувязанно с задачей определения частотных характеристик и напряженно-деформированного состояния при переменной во времени внешней нагрузке.
Проведенное усложнение математической модели системы управления несущим винтом, отражающее податливость упругого торсиона воздействиям системы переменных нагрузок, значительно расширит возможности конструктора при проектировании несущей системы.
Созданная методика численного моделирования позволяет на ранней стадии проектирования оценить взаимодействие жестких и податливых звеньев механической системы под влиянием заданных нагрузок, разработать методики определения частотных характеристик системы и эксплуатационный диапазон по режимам нагрузок. При необходимости конструктор имеет возможность моделирования с определенной степенью точности некоторых режимов стендовых и летных испытаний до построения натурных стендов и летных испытательных образцов.
ЛИТЕРАТУРА
1. Бате К., Вилсон Е. Численные методы анализа и метод конечных элементов / Пер. с англ. А. С. Алексеева и др.; Под ред. А. Ф. Смирнова. — М.: Стройиздат, 1982, 448 с. ил. — Перевод изд.: Numerical methods in finite element analysis / K.—J. Bathe, E. L. Wilson (1976).
2. Белоцерковский С. М. Тонкая несущая поверхность в дозвуковом потоке газа. — М.: Наука, 1965, 244 с.
3. Белоцерковский С. М., Скрипач Б. К., Табачников В. Г. Крыло в нестационарном потоке газа. — М.: Наука, 1971, 768 с.
4. Братухин И. П. Проектирование и конструкция вертолетов. — М.: Оборонгиз, 1955, 360 с.
5. Втулки несушцх винтов вертолетов. Обзор № 393 ЦАГИ им. Н. Е. Жуковского. Составитель Лернер М. А. // ОНТИ ЦАГИ, 1972, 63 с.
6. Голованов А. И., МитряйкинВ. И., Шувалов В. А. Исследование напряженно-деформированного состояния торсиона бесшарнирного несущего винта в геометрически нелинейной постановке // Вестник МАИ. Т. 15. № 5, стр. 44 — 53.
7. ЗенкевичО. Метод конечных элементов в технике. — М.: Мир, 1975, 270 с.
8. Миль М. Л., Некрасов А. В., Браверман А. С., Гродко Л. Н., Лей-кан д М. А. Вертолеты. — М.: Машиностроение, 1966. Кн. 1. — 455 с., 1976. Кн. 2. — 424 с.
9. Морозов В. И. Математические модели динамики аэроупругого летательного аппарата. Исследование авиационной техники с помощью ЭВМ. — В кн.: Труды ВВИА им. Н. Е. Жуковского, 1981, вып. 1310.
10. Образцов И. Ф., Савельев Л. М., Хазанов Х. С. Метод конечных элементов в задачах строительной механики ЛА. — М.: Высшая школа, 1985, 392 с.
11. C г a i g R. R. and Bampton M. C. C. Coupling of substructures for dynamics analyses // AIAA J. 1968. 6(7):1313—1319.
Рукопись поступила 7/V 2009 г.
