CC BY
41
УЧЕНЫЕ ЗАПИСКИ Ц АГ И Том XVIII 1987
№ 2
УДК 518: 512.35
ПРОБЛЕМА СОБСТВЕННЫХ ЗНАЧЕНИЙ МАТРИЦ В ЗАДАЧЕ О ФЛАТТЕРЕ
Новый метод решения полной проблемы собственных значений неконсервативной системы, учитывающий особенности уравнений флаттера и ориентированный на решение задачи о флаттере. Применяется дифференциальная интерполяция по параметру, каким в задаче о флаттере является скорость потока.
Алгоритм [1], рассчитанный на решение задач о флаттере с г=30 степенями свободы, как показала десятилетняя расчетная практика, устойчиво работает только до г=15. Поэтому расчеты на флаттер, сводящиеся в своей алгебраической части к проблеме собственных значений матрицы 2/--ГО порядка, где г — число степеней свободы, следует проводить, пользуясь не прямыми методами, а косвенными, учитывающими особенность решаемой задачи. Достаточно отметить, . что для уравнения флаттера:
Матрица А почти наполовину пустая, а комплексные частоты 5 = 6-И'о> слабо отличаются от чисто мнимых: |6|<С(», т. е. система имеет ярко выраженные резонансные свойства.
Целью предлагаемой методики является построение траекторий корней в (годографов частот) при изменении параметра V от нуля до заданного максимального значения. На это обстоятельство обращается внимание в книге Уилкинсона ([2], гл. 9, § 61) как на ключ для создания методики; а именно, решение для каждого нового значения параметра V ищется как уточнение решения, полученного на предыдущем шаге по V.
В предлагаемом методе за основу берется варьирование параметра и (поскольку этот параметр обычно выполняет роль скорости полета, то далее мы будем его называть просто скоростью).
В. Г. Буньков
ои + си + V' ви + V йи = 0, где С, С, В, £> — матрицы г-го порядка, имеем:
эа
Особенности предлагаемого метода:
— на каждом шаге по скорости за исходное берется не просто предыдущее решение, а поправленное с учетом производной, т. е. применяется дифференциальная интерполяция;
— удается остаться в рамках матричного уравнения г-го порядка, не переходя к матрице 2г-го порядка;
— разработан специальный алгоритм для преодоления затруднений в случае бифуркаций (встреча двух траекторий в точке кратности комплексных корней);
— в качестве исходного значения скорости берется и = 0, соответствующее колебаниям в пустоте.
Задача о колебаниях в пустоте приводит к уравнению п-то порядка для консервативной системы:
вХ + СХ = 0 или 6^ = X С*; X = оу* , (1)
где в, С — симметричные матрицы: жесткостная и инерционная.
Для любого вектора выполняется условие:
0; гтсг>о.
Способ определения корней (частот) уравнений (1) хорошо известен и сводится к следующей схеме:
а) С = 14.; I—треугольная; =
б)АУ=ХУ\ X = МУ; Л = УИтО/И — симметричная;
в)£) = (2тЛ<3 — трехдиагональная (метод Гивенса); DZ = XZ; Л' = М(}г;
а1 р2 0 0 0
Ра «2 Рв 0 0
0 Рз а3 р4 0
г) собственные числа трехдиагональной матрицы вычисляются методом бисекций (последовательность Штурма) с подсчетом числа перемен знака.
Определение собственных векторов представляет более трудную задачу, чем определение корней. Для ее решения предлагается вариант метода обратной итерации [2], который можно назвать методом управляемой обратной итерации. Анализируется расположение корней в порядке неубывания Л1<А,2< ... Хг\ г*сп. Каждый корень %г возмущается на величину ¿¿ = 8^ (е1 — 0,01), где А» — расстояние до ближайшего корня, не считая кратных (критерий кратности: \Хг—Я,г+11 <езЯ-г, б2~Ю~4) и вместо уравнения (/)—Х{Е)1 = 0 решается уравнение = 0;
= Д—ХЕ, где X—Хг 4“ бг — возмущенный корень. Последнее уравнение решается с помощью итерационного процесса
0\21+х=г], у=0, 1, 2, ... , (2)
где в качество исходной итерации берется набор псевдослучайных чисел, например, Мь = 20 + к+ 10- (—1,05)*; Zo= («1... ип).
Уравнение (2) удобно решать в два приема, в первом из которых трехдиагональная матрица приводится к двухдиагональной:
рг = ь ; /1 т, о о ... \
0 1 ъ о . . . а = (ах... а„)т; Ь = (Ь1 ... ЬпУ; \...................../.
Алгоритм такого преобразования можно записать так:
1. С1 = а,- К; Та = Рз/^1 ; с = \с11.
2. Цикл к, — 2, ... , п: ск = а.к — X - ^-1 с: = с+ |сА|; = Р*+1/сЛ.
3. Нормировка: ск: = ск1с\ к— 1, ... , п.
4. Решение системы: = а^с^, Ьк — (ай — Ьк-\!с)!ск; к = 2, ..., я.
5. = г„ = Ь„ — 1к^к+1\ к — п- 1, ... , 1.
6. Нормировка: /=^|г*|; ■£* : = £*//, к — 1,... , п.
7. Критерий сходимости^ |гй —ай |<е3 или ^1+ а*К£з-
Для машинной точности е = 2х10"12 рекомендуется ез=Ю~10. Найденный таким образом вектор Ъ для %г ортогонализируется ко всем предыдущим ги...,
чг-я*г-, я = \г,... г^]-, г,=г-яш.
После этого вектор Zi нормируется: ZJ Еь=\. Для получения каждого нового вектора следует обновлять исходную итерацию Z0, например, циклической перестановкой компонентов прежнего набора случайных чисел.
Таким образом, решается задача о колебаниях системы в пустоте. Для системы порядка п=20-ь85 рассчитывается г — 5-н20 тонов колебаний, которые берутся для расчета на флаттер, т. е. делается переход от системы с матрицами «-го порядка: Ап = Оп, Сп, Вп, Оп к г-му через базис Б [г+п], состоящий из собственных векторов:
А = БАпБт; Бт = (Х1г..., Хг) .
Если векторы нормированы по кинетической энергии: БСпБт = Е, то уравнение флаттера г-го порядка приобретает канонический вид:
в, и + Д 0 + 0 = 0; G1^G + v*B; А=Д> + «Я» (3)
здесь для общности добавлена матрица трения О0. Матрица жесткости
2 2
в является диагональной, состоящей из квадратов частот: а)1<ш2<
< ... .
Приступаем к основной задаче: расчет колебаний неконсервативной системы, описываемой уравнением (3).
Будем считать, что нам известно решение и = иедля скорости V, т. е. даны комплексные частоты 8и = &и + Шк и формы ик=ик + Шк колебаний ,(к= 1,..., г) и требуется найти решение для новой скорости v + v/, где V' — шаг по скорости. Система (3), вообще говоря, имеет 2г корней, т. е., кроме указанных г корней, имеются еще г комплексно-сопряженных: Ьк-~ ¿<ок, ик — Юк. Однако возможен случай, когда
вместо одной или нескольких пар комплексно-сопряженных корней имеют место столько же вещественных — в таком случае придется рассчитывать не г траекторий, а несколько больше, но сущность описываемого ниже метода от этого не изменяется, так как каждая траектория
рассчитывается отдельно.
Итак, имеем частоту 5 и вектор и. Кроме того, нам понадобится еще и вектор V транспонированной системы:
О! У + Ц]У+ У = 0 . (4)
Уравнения (3) и (4) адэкватны следующим уравнениям в форме
Коши:
А0 = зи ; Ат У=5У,
Воспользуемся известным соотношением для определения производной частоты 5 [3]:
^ т ¿^4 —
Ут — и
сш
Ут V '
( у _д_(-Ог-вЛ (з1/\
\ — бт V ) ди 1 Е 0 / V и )
— б1- V
И
«у
и
{зуТ-.-У1 Од
-И
о
-X!/
о
х-1
в
■аВ) и
(з Ут; -Кт О
■>(«)
Vт (й1 — Е)1/
(5)
а = у. г»
|Х—1
Для исходного варианта по скорости (1^=0) частоты и формы £-го тона известны: в = гсо^; £/— У= (0... 1... 0)—единица на к-м месте. Поэтому первый шаг для частоты 5 при подстановке в (5) будет:
Ькк
2«б
(6)
где Ькк, Лкк — диагональные элементы матрицы В, Б. Что касается матрицы трения £>0, то, во-первых, трение предполагается малым с декрементом 6<0,1, во-вторых, в исходном варианте следует его'учесть, взяв вместо чисто мнимой частоты Ши поправленную на трение: 5 =
— — -|- I о)А .
Для нулевых и очень малых частот, характеризуемых неравенст-
вом: <вй < е3 (оср; й = 1 ... г0; (ос =
’ шк
формула (6) не годится и
заменяется точным расчетом для подсистемы, состоящей из г0 тонов. Вопрос о нулевых тонах в потоке может быть решен отдельно.
Приращение частоты Дя для произвольной скорости V согласно (5) определяется по такой схеме:
£/, = «¿7; и2 — «£/х ; £/3 = б, £/-£/,
ил = 01и + 01и1 + и2; иь = ю'Ои2 + о.'Ви1; а'-®'*®1
А8 = Угиб/Ут и»-, о == | £/4 |/о
ср
Здесь также написано выражение для относительной ошибки а, характеризуемой невязкой при подстановке решения в исходное уравнение (3). При этом полагается, что вектор нормирован: | С/1 = 1 (на практике удобнее другая нормировка: наибольший по модулю компонент вектора и равен единице).
Приращение частоты Ая позволяет в линейном приближении получить ожидаемую частоту в0. Частота 5 и векторы и, V нового варианта по скорости V являются уточнениями предсказанной частоты 50 и векторов прежнего варианта. Определяется поправка I к частоте:
Из уравнения (3) получаем:
или
Е ¿7= О
} (7)
где
6* — (@1 Ч- 5о 4* *о2 Е) ! В* — -{- 2х0 Е .
Уравнение (7) можно переписать так:
(Ї.Е-Н) и = У?и; АУ*=\У\ Л = ^ ; У =
Я = С;1 ; /=• = //£>*.
Искомое число Я является наибольшим собственным числом матрицы А, превосходящим по модулю остальные на много порядков и поэтому легко определяемым итерационным (степенным) процессом. Чтобы избежать осложнения при слишком хорошем совпадении 50 и 5, можно ввести возмущение ожидаемой частоты 50:
где — ближайшая по модулю ожидаемая частота, не считаякратных: | Я; — 50 К £2 “ср- Вместо итерации матрицей А 2г-го порядка предлагается итерационный процесс с матрицей г-го порядка:
где хт—максимальный по модулю компонент вектора I (т — номер компонента); —ит/хт — отношение т-х компонентов векторов ¿7,-1 и /. В исходной итерации (у=1) задается £о=0, а в качестве вектфа и0 берется собственный вектор из предыдущего варианта по и. Условие сходимости итерационного процесса (8) записывается так:
Искомая комплексная частота в случае сходимости (9) или превышения максимального разрешаемого числа итераций Vmax, а также вектор ІІ определяются так: £ = 50—и =1)ч .
Параллельно с процессом (8) производится процесс для вектора V:
Так, идя с шагом V' по скорости и, строятся годографы комплексных частот 5 до тех пор, пока не возникнет случай кратных корней (близких). Такой случай, называемый бифуркацией, характеризуется встречей двух траекторий, идущих навстречу друг другу. После встречи в точке У* при некотором значении скорости v — v¡¡: траектории резко расходятся в противоположных направлениях, но уже под прямым углом к прежним (рис. 1). В окрестности бифуркации поведение двух корней можёт быть описано формулой:
| |/соср < £3 ; V > 2 .
(9)
J =ЯТ [£[ У,_! + (2я0 - 1Л_,] ; V, = Г\хт .
«1,2=5*+ сУ v — v9+o(v — vй;).
7 «Ученые записки» № 2
97
Расстояние между корнями Л = 1—^21 в этом случае убывает по параболическому закону: к — | 2с ]/г>—|, а последний (к-й) шаг по скорости V перед бифуркацией характеризуется неравенством:
Ь.\-\ — 2И1 > 0 . (10)
Этот признак служит для замены дифференциального прогноза бифуркационным:
«1,2 (к + 1) = 2/к —/к-1 + У2§*к — ¿к-\ ,
/и = [^1 (£) + 52 (£) ]/2 ; gk = [52 (к) — {к) ]/2 .
Следующий шаг после бифуркации может оказаться тоже удовлетворяющим неравенству (10), но на этот раз расстояние /г растет, а не убывает. Поэтому к критерию (10) следует добавить критерий убывания расстояния:
Бифуркация предсказывается формальным перебором всех пар найденных на текущем шаге корней. Векторы для итерационного процесса при бифуркационном прогнозе берутся в любом сочетании из имеющихся двух векторов до бифуркации.
Особым случаем бифуркации является превращение комплексного корня в пару вещественных (рис. 2). При этом возникает дополнительная траектория частоты. Такой случай характеризуется неравенством:
¿1 = <4_1 — 2о)*> 0 , а прогнозируемые вещественные корни определяются так:
$1, 2 = 8* +V (I •
Описываемый процесс бифуркации пригоден и в случае ситуации, близкой к бифуркации (пунктирные кривые на рис. 1), тем более, что точная бифуркация комплексных корней на практике встречается редко.
Предлагаемый метод не встречает затруднений в случае кратных частот в начале траекторий корней (в пустоте), так как все векторы ортогональны.
В результате расчетов установлено:
Расчет г—24 тонам колебаний в пустоте для системы п = 60-го порядка занимает на ЭВМ 150 с с ошибками 0=10“® для низших тонов и а = 10_3 — для высших.
Расчет комплексных частот в потоке с г = 24 занимает 25 с на каждый вариант по скорости. Несмотря на то, что получающиеся при этом ошибки очень малы: 0=1О~10, порядок г = 24 является, по-видимому, предельным, выше которого начинается численная неустойчивость результатов, проявляющаяся в том, что, несмотря на ничтожную невязку о=10-10, ошибки в комплексных частотах 5 = 6-И'со, особенно в б, достигают нескольких процентов, что может привести к неправильным выводам о критической скорости флаттера.
Обычное число шагов 10-т-30 при вариации скорости от нуля до максимальной всегда приводило к устойчивому расчету всех ветвей годографов комплексных частот. Дробление шага еще втрое не нарушало сходимости. Значит, рекомендуемый шаг по скорости составляет 0,01—0,10 от максимальной.
ЛИТЕРАТУРА
1. Буньков В. Г. Проблема собственных значений матриц в расчетах на флаттер. — Ученые записки ЦАГИ* 1975, т. IV, № 2.
2. Уилкинсон Дж. X. Алгебраическая проблема собственных значений. — М.: Наука, 1970.
3. Ф а д д е е в Д. К., Фаддеева В. Н. Вычислительные методы линейной алгебры. — М.: Физматгиз, 1963.
Рукопись поступила 9/Х 1985 г.
