CC BY
66
УЧЕНЫЕ ЗАПИСКИ Ц А Г И
Т о м XII 19 8 1 №2
УДК 629.735.33.015.4:533.6.013.422
ПРИМЕНЕНИЕ МЕТОДА Д-РАЗБИЕНИЯ ДЛЯ ПОСТРОЕНИЯ АЛГОРИТМА РАСЧЕТА НА ЭВМ УСТОЙЧИВОСТИ ЛИНЕЙНЫХ СИСТЕМ И СИСТЕМ С ЛОКАЛЬНЫМИ НЕЛИНЕЙНОСТЯМИ
А. В. Крапивко
Применительно к задачам флаттера рулевых поверхностей и шимми колес самолета рассмотрен алгоритм расчета периодических решений для линейных систем и систем с локальными нелинейностями, основанный на методе Д-разбиения и частотных характеристиках динамической жесткости. Приводится пример расчета устойчивости движения ориентирующихся колес стойки шасси самолета, оснащенной демпфером шимми с линейной и нелинейной характеристиками силы сопротивления.
При изучении явлений флаттера рулевых поверхностей и шимми колес шасси самолетов необходимо учитывать ряд существенных особенностей реальных конструкций, которые в значительной степени ограничивают применение широко известных методов расчета [1 — 3] и требуют их совершенствования и развития. Одной из этих особенностей является большое разнообразие конструктивных схем систем управления, что затрудняет их единообразное математическое моделирование и, как правило, приводит к значительным изменениям алгоритма решения поставленной задачи при переходе от одной конструкции к другой.
Вторая существенная особенность исследования устойчивости движения реальных конструкций заключается в необходимости учета ряда нелинейных факторов, среди которых наиболее важными являются нелинейности характеристик механизмов управления (бустеров, электрогидравлических следящих приводов), противофлаттерных демпферов или демпферов шимми, люфты и силы сухого трения между подвижными частями конструкции, а также нелинейности аэродинамических сил или сил взаимодействия пневматика с поверхностью земли.
В настоящей статье рассматривается алгоритм расчета на ЭВМ границ области устойчивости для линейных систем и параметров предельных циклов колебаний для систем с локальными нелинейностями, который разработан на основе метода Д-разбиения в плоскости одного комплексного параметра [1]. Благодаря единой формализации вычислительных процедур для систем с различным числом степеней свободы алгоритм позволяет в значительной мере преодолеть указанные выше трудности. Основная особенность алгоритма состоит в использовании для Д-разбиения плоскости комплексной величины приращения к элементу ^.-матрицы, соответствующей исходной системе дифференциальных уравнений.
Анализ устойчивости нелинейных систем проводится на основе метода гармонической линеаризации.
1. Определение границ области устойчивости для линейных систем. Рассмотрим задачу об определении границ области устойчивости системы, математическая модель которой может быть представлена в виде системы линейных
уравнений с постоянными коэффициентами. С помощью преобразования Лапласа приведем эту систему уравнений к системе однородных алгебраических уравнений с комплексной матрицей размером N X N
К (1)0 = 0 (1)
с элементами /Су* (А). Здесь X = 8 + г'ю— комплексная частота колебаний (3—показатель колебаний, и — частота колебаний); (2 = (д 10........9л?о)т — вектор ком-
плексных амплитуд (форма колебаний). Верхний индекс „т“ означает транспонирование.
Введем в рассматриваемую систему дополнительный элемент, изменение параметров которого эквивалентно приращению Д/Су* (А) к элементу /Су* (X) матрицы К(X). Для построения границ области устойчивости на плоскости параметров этого элемента выберем комплексную величину Д/С/* (X) в качестве плоскости для Д-разбиения. Связь этой величины с элементами матрицы /С(Х) находится из условия существования нетривиальных решений уравнения (I):
где с!е1 'К (X) — определитель матрицы К (X) при Куй = 0; (Ш А]ц (X) — алгебраическое дополнение элемента (X) матрицы АГ|(Х); Н0 (X) и (X) — характеристические полиномы соответственно исходной математической модели и этой же модели, вектор координат которой не содержит ?£о'ГО компонента.
Построив годограф при значениях 0^ш<;оо, получим ото-
бражение мнимой полуоси <» ;> 0 плоскости корней характеристического уравнения системы на плоскость приращения Д/Су*- При исследованиях флаттера рулевых поверхностей и шимми колес из физических соображений, как правило, может быть определена степень устойчивости (или неустойчивости) системы при нулевых или бесконечно больших значениях ДАГ;й(ги)). которыми определяется, например, динамическая жесткость системы управления. Это обстоятельство существенно упрощает выделение области устойчивости на плоскости Д/С;*.
Обозначим через о>кр значение со, при котором 1?е Д/(д (/<окр) = 0 или 1т Д/Су* (гмКр) = 0. Соответственно обозначим
(1т Д/Су* (г’“кр) при Ие Д/Сд (гмкр) = 0,
/Кр)= (3)
I Яе Д/Су* (г«кр) при 1т Д/Сд (гшКр) = 0.
Если исследуемый элемент системы такой, что соответствующая ему величина ДК;ь(1и>) имеет только действительную или только мнимую часть, то его критические с точки зрения устойчивости значения могут быть определены с помощью функции /(<йкр).
Пусть известно, что при Д/Суй (Ло) = 0 система устойчива. Тогда значение функции /(сокр) при пересечении годографом Д/С;-Й(гш) мнимой или действительной оси с возрастанием его аргумента при увеличении а> определяет максимально допустимое (по модулю) для устойчивости значение ДК/йО'мкр) и> наоборот, если известно, что при | ДЯуй (/<») | = оо система устойчива, то значение функции /(шкр) при пересечении годографом ДАТ;* («“) действительной или мнимой оси с убыванием аргумента при увеличении ш определяет минимально необходимое (по модулю) для устойчивости системы значение ДД’уд; (гикр) [1].
Варьируя параметр задачи рт, можно также определить его критические значения в зависимости от параметров дополнительного элемента ДКу*.
Из соотношений (2) и (3) видно, что эффективность расчета функции /(о>кр) определяется процедурами вычисления годографа Ди корней уравнений Ие ДАСу* (гш) = 0 или 1ш Д/Су* (га>) = 0 из (3). Из (2) следуют два возможных пути расчета годографа Д/С/* (гм): первый путь основывается на вычислении определителей матриц 'АГ(Х) и Лу* (X), второй — на вычислении характеристических полиномов #0 (X) и Н1 (X). Из них более рациональным является первый путь,
так как он не требует формирования сложных матриц порядка необ-
ходимых ДЛЯ построения ПОЛИНОМОВ #о(Х) И Нг (X), и позволяет без изменения структуры матрицы \К (X) эффективно использовать в расчетах различные передаточные функции для формирования любого ее элемента.
Для вычисления определителей с1е1 'К (X) и с!еМд(Х), а также форм колебаний системы на границе устойчивости в разработанном алгоритме используется декомпозиция (¿¿/-разложение) соответствующих несингулярных матриц [4]. Матрицы 'КО') и А]к (X) представляются в виде произведения нижней треугольной матрицы /. (X) с единицами по главной диагонали и верхней треугольной матрицы и (К) : 'К (X) = /. (X) ¿/ (X).* Из этого равенства получаются простые рекуррентные соотношения для вычисления элементов матриц ¿(X) и и (К).
Очевидно, что det 'К (X) = det L (X) det ¿7(X) =
Uu (X), где UuQ.) — диагональ-
i=i
ные элементы матрицы <У(Х).
Для обеспечения сходимости процедуры разложения матрицы ее следует предварительно балансировать путем приведения исходной системы уравнений к безразмерному виду, масштабировать с помощью масштабного множителя &о К (X) = &0 К (X), а затем, следуя [4], использовать прием, называемый частичным выбором главного элемента.
Для вычисления корней шКр уравнений Ие Д/Су* (/со) = 0 или 1т ДЯу* (/со) = 0 из (3) используется метод секущих [5], который для однозначных и непрерывных функций НеДАГуй(гш) и 1т Д/С,* (/ш) имеет сверхлинейную скорость сходимости. Итерационный процесс метода секущих заканчивается на (к + 1)-й итерации при выполнении двух условий:
£ol/i(“b+i)IO
f(°>k+1)+/(“*)
< s2>
где /, (сол+1> = Ие Д/Су*(го>А+1) или Л (сой+1) = 1т Д/СуА (/соА+1); и е2 —заданные величины погрешностей.
В качестве искомых значений корня сокр и функции /(сокр) принимаются следующие величины:
°Кр — 0.5 (°>£ + 1 + “ft)! / (МКр) ■— [/(<0А+1) + /(“ft)].
(4)
Если корень X = /сокр является простым (некратным) собственным значением, то в качестве базисного минора матрицы К (X) можно выбрать регулярную подматрицу Л у* (X), разложение которой получается при вычислении приращения Д/Суй (/со). Заменим уравнение (1) эквивалентным уравнением:
= 0,
Кы-1. 1 (г'“кр); ^,V—1, 1 (г"кр) ’ Qn- I
^i, jv—1 (г"“кр); ^1. 1 (гшкР) Qi
где
ATjy__X, /V—! (2шкр) — Лyft (/сокр); (гсокр) — ACyft (/“кр) ^Kjk (*“кр); Qi — Ч\
&Д1_1 1 (г’“кр). (/“кр) — внебазисные й-й столбец и /-я строка матрицы
К (/с»кр)-
Тогда вектор 9 = (9Л/—1> *?)т (форма колебаний) определяется из решения следующей системы неоднородных алгебраических уравнений:
К
N—1, Л7-1
(гшкр) Qn—1 — &N—1, 1 (г“кр) Ql-
(5)
Но так как ___1, л'’_1 (г“кр)== Лу* (¿соКр)—I*(/с°кр) £/дг_1,/V—1 (г(йкр)> то
решение уравнения (5) сводится к последовательному решению двух систем уравнений с треугольными матрицами:
¿w-l. JV—1 (г“кр) ^N—1 — ^N— 1, 1 (г“кр) <7>
Uh-1, N—1 (г“кр) Яц-1 = 1 •
(6)
Для окончательного суждения о точности вычислений шкр и ф использует ся следующая интегральная оценка а:
II К (/»*р) <? II
(7)
“кр II 0 II
где квадрат нормы вектора равен сумме квадратов модулей его компонентов.
2. Определение параметров предельных циклов! колебаний для систем с локальными нелинейностями и оценка их устойчивости. Применим изложенный алгоритм для приближенного исследования автоколебаний автономных систем с локальными нелинейностями, т. е. для определения периодических режимов колебаний (предельных циклов) и оценки их устойчивости. Будем рассматривать только такой класс нелинейных систем, когда нелинейности образуют не более чем две группы, каждая из которых определяется одной из обобщенных координат.
Для решения задачи параметрического исследования автоколебаний типа флаттера органов управления и шимми колес самолета воспользуемся методом гармонической линеаризации. Будем предполагать, что применение этого метода для указанных задач обосновано [6].
Комплексная матрицы К (г'м, нелинейной системы однородных уравнений будет содержать два элемента К"к (««>) + (г®, <7Й 0) и А^г(г'<о) +
+ (¿м, <7г0), зависящих от компонентов дк0 и <7; 0 , вектора решения 0 (гм) и
частоты м (индексы „л“ и ,н‘ обозначают линейную и нелинейную части элемента матрицы). Для ряда важных случаев нелинейные части этих элементов могут быть получены в виде аналитических зависимостей, а в более сложных случаях вычислены путем гармонического анализа численного решения соответствующей системы нелинейных уравнений [7]. Для реальных объектов величины (гш, ок0) и /('"¿(г'м, q[0) могут быть определены по результатам частотных испытаний [8]. В последних двух случаях величины (г’“> <7Й 0), /С”г(г'м, <7го) задаются для дискретного ряда значений двух аргументов м, о или д[(). Для вычисления значений этих величин в промежуточных точках удобно использовать сплайн-функции [9].
Параметры предельных циклов колебаний системы могут определяться по следующему алгоритму:
1) для фиксированных значений аргументов м и <7;о вычисляются значения
величины (гм, qlй) по соответствующим формулам или с помощью сплайн-
функции;
2) с помощью соотношения (2) определяется значение частотной характеристики (гм, <7го);
3) из уравнения (6) находится вектор <? (гм) и, в частности, модуль его к-то компонента.
Таким образом, для различных значений амплитуд 9;(), 0 и частоты « мо-
жет быть вычислено семейство частотных характеристик Д/<уА(г'м, д[0) и модуль цкй, необходимые для существования предельных циклов колебаний. Варьируя параметр рт в равенстве Д/Су* (гм, д[0) = (гм, ^ 0), можно найти искомые характеристики предельного цикла колебаний в зависимости от рт.
В общем случае, когда действительная и мнимая части /Су* (гм, Якй) зависят от двух аргументов со и ^0, определение параметров предельных циклов колебаний из указанного выше равенства является весьма трудоемкой задачей. Эта задача эффективно решается в двух следующих часто встречающихся на практике случаях: коэффициент жесткости а^(дк0) или демпфирования (^^ 0) системы управления зависит от амплитуды колебаний <7^0 органа управления (люфт в проводке управления или демпфер с нелинейной характеристикой силы сопротивления).
В этих случаях
и соотношения (2) — (6) разработанного алгоритма позволяют вычислить необходимые для существования предельных циклов колебаний значения функции /(“кр><7й0> <7;о) и найти точки ее пересечения на плоскости с заданными функциями а.к(,ц0) или к.к(Чк0).
Погрешность расчета параметров предельных циклов колебаний определяется, как и в линейном случае, интегральной оценкой с (7).
Устойчивость предельных циклов колебаний можно приближенно оценить, сопоставляя расчетные значения величины Д(г'ш, 0, о) и заданные ее зна-
чения /("¿(гм, <7Й0). Для двух указанных выше случаев условия устойчивости могут быть сформулированы следующим образом: предельный цикл будет устойчивым, если для соответствующей линеаризованной системы значения функции /(“кр. 0. ?го) являются минимально необходимыми для устойчивости значения-
ми и при увеличении амплитуды дк0 от ее критической величины эти значения остаются меньше заданных величии (г®, <?А0); если значения функции/(“Кр. ^/го’ <?, о) являются максимально допустимыми для устойчивости линеаризованной системы значениями, то для устойчивости предельного цикла необходимо, чтобы с ростом амплитуды 0 от ее критической величины эти значения были больше заданных величин /Су* (гш, Чьо)-
3. Пример расчета. Проиллюстрируем применение разработанного алгоритма на примере исследования устойчивости качения ориентирующихся колес стойки шасси самолета, оснащенной гидравлическим демпфером для предотвращения шимми. Для этой цели воспользуемся полученной в работе [10] системой
линейных уравнений шимми и для упрощения пренебрежем в них соотношениями, связанными с блокировкой колес. Выберем вектор обобщенных координат в виде ф = (ф, 0, Х/г, <э)т, где ф и 0 — углы поворота стойки шасси относительно осей Ох и Оу (рис. 1,д); X и <р—смещение центра контакта шины Ісолес в боковом направлении и угол поворота ее зоны контакта. Тогда соответствующая матрица Я’(г'“) в безразмерной форме имеет вид:
- - Тх ш3 +“Сф 72 + Т ¿2; _1ХУ»2 + ій 1 + т, 0;
— .— —/дгу — гУ1К] —/„ м2 -¡- £> (г'ш) -(- гш/г22; і; Ь\
К (го>) = _ ___ __
/о>(1 + /); V + Ш\ ¿'о>; К;
_ — ^ К; г'ш; —а К; р"К -Ь -
где 1Х = ТХ птпг-, 1у — ТуПтг2, 1ху = /ху л/гаг2 — моменты инерции установки шасси относительно соответствующих осей; /* = /* /иг2 — момент инерции колеса относительно оси его вращения; V = —скорость движения самолета; к =
= ка, Ь — Ьаг2 — радиальная и пяточная жесткости шины; с^ = с^па — боковая жесткость установки шасси; а = аг2, ¡3 = Зг, у = уг — кинематические характеристики шины; £ = їг — вынос колеса; / = ¿г — приведенная длина стойки, и = = ш]Аа//?г; а, т, г и п— боковая жесткость шины, масса, радиус колеса и число колес (масштабные коэффиценты).
Второе слагаемое элемента АГзгО“) является частотной характеристикой динамической жесткости связи между ориентирующейся (1) и неподвижной (2)
где
амплитуда угла 6.
_ _ Мв (г<о)
частями стойки шасси (рис. 1, б): В (гю) — паГ2$.~
Представим зависимость между скоростью относительного перемещения поршня
Мо ‘
демпфера и усилием его сопротивления Мд [10] в виде М$ =/г ^0 — ) > ГДе
11 = Тмг2~угат — коэффициент сопротивления идеального демпфера, с9 = лаг2 — жесткость стойки на кручение. Из этого уравнения при 0(£) = 0оегт* и Мд (() = = ./М6 (гм) имеем
Re D (гм) = .
, Im D (г'ш) =
1 +
= /гм где (1) = - .
с0
По соотношениям (2) — (4) были вычислены значения функции
/(“кр. /г) = ~ 1ш Д/С22 (г«Кр) = -Л22 мкр >
которые позволяют оценить влияние на шимми дополнительного демпфирования й22, возникающего главным образом из-за момента Мтр сухого трения между подвижными частями стойки при ее крутильных колебаниях (см. рис. 1, б).
Представленные на рис. 2 зависимости тс/г22юкр от Л показывают, что установка демпфера обеспечивает отсутствие шимми, если коэффициент сопротивления демпфера Л находится в интервале значений, при которых функция /(“кр> А)<0. Наличие в стойке дополнительного демпфирования повышает запас устойчивости по Ъ.
В случае квадратичной зависимости силы сопротивления демпфера от скорости относительного перемещения его поршня
Ma = Лк 0 —-
Мд
Cft
sign 0
Щ
Со
где hK — h^nmr2 — коэффициент сопротивления нелинейного демпфера.
С помощью гармонической линеаризации можно получить следующие формулы для вычисления составляющих динамической жесткости:
Re D (iu>, 0О) = £)д; Ira D (iu>,
в.- Ÿ
Зл
i6 ;
\2 =2 )
i = A,Vl -Z>| ; Зтс
16 ;
= hv w ,—
где Ci) - — Y0O.
V CB
Для данного случая зависимость функции /(«кр. ^о. Лк) = "кр от ПРИ V = 1,31 показана на рис. 3 сплошными линиями. Там же представлена соответствующая зависимость шк.р = ыкр \гт1а . Если величинами дополнительных сил трения в стойке шасси пренебречь, то точки пересечения функции / (60) с осью абсцисс определяют частоты и амплитуды искомых предельных циклов колебаний (рис. 4, Мтр = 0).
В случае значительных по величине сил трения в стойке шасси их влияние на устойчивость можно приближенно оценить, рассматривая зависимость момента сил сухого трения Мтр от скорости 0 в виде идеальной релейной характе-
4 Л?тр
ристики. Тогда ТОЧКИ пересечения функции ЛЫЙ22КВ (0О) = 0-— (рис. 3, штриховые линии), где Л|*в (0О)—коэффициент гармонической линеаризации релейной характеристики, с функцией ла>кр Л22 (0О) определяют частоты и амплитуды предельных циклов колебаний с учетом сил сухого трения (рис. 4, Л4тр ф 0).
----гщмотеш
пинеаризшия \
15
Рис. 4
Из сопоставления представленных на рис. 2 и 4 результатов следует, что локальные нелинейности тийа сухого и квадратичного трения существенно влияют на устойчивость движения ориентирующихся колес самолета. Проведенное исследование с учетом этих нелинейностей указывает на существование устойчивых и неустойчивых циклов колебаний. Сравнение результатов расчета по предложенному алгоритму с соответствующими результатами численного интегрирования (рис. 4) указывает на удовлетворительную для инженерных расчетов точность метода гармонической линеаризации при исследовании шимми с учетом рассмотренных нелинейных зависимостей.
В заключение отметим, что время расчета одной точки границы области устойчивости для линейной модели шимми с помощью разработанного алгоритма сокращается не менее чем в пять раз по сравнению с временем решения этой же задачи с помощью алгоритма, использующего схему Рауса и решение полной проблемы собственных значений матрицы [2, 3].
ЛИТЕРАТУРА
1. Ней мар к Ю. М. Устойчивость линеаризованных систем. ЛКВВИА. 1949.
2. Буньков В. Г. Полная проблема собственных значений матриц в расчетах на флаттер. „Ученые записки ЦАГИ\ т. VI, № 2, 1975.
3. Соболев Е. И. Применение цифровых вычислительных машин для определения критической скорости флаттера систем со многими степенями свободы. Труды ЦАГИ, вып. 949, 1968.
4. Уилкинсон Дж., Р а й н ш К. С. Справочник алгоритмов на языке АЛГОЛ. Линейная алгебра. М., .Машиностроение“, 1976.
5. О с т р о в с к и й А. М. Решение уравнений и систем уравнений. М., Изд. иностр. лит., 1963.
6. Попов Е. П. Прикладная теория процессов управления в нелинейных системах. М., „Наука“, 1973.
7. К р а п и в к о А. В. Исследование эффективности демпферов для гашения колебаний органов управления. Труды ЦАРИ. вып. 1871, 1977.
8. Forschung H. Einfluß servomechanischer Steuerungs und Stabi-iitätssystem auf das Flatterverhalten von Flugzeugen. Z. Flugwiss. 21 (1973), Heft 1.
9. H a r d e r R. L., Desmarias R. N. Interpolation using surface splines. „J. Aircraft“, II, vol. 9, N 2, 1972.
10. Г о з д e к В. С. Устойчивость качения сблокированных ориентирующихся колес самолета. Труды ЦАГИ, вып. 1196, 1970.
Рукопись поступила 25/VI 1979 Переработанный вариант поступил 2411II 1980
