Применение нового математического аппарата «Одномерные спектральные портреты матриц» к решению проблемы аэроупругих колебаний решеток лопастей Текст научной статьи по специальности «Физика»

Том X Ь

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

№ 6

УДК 518:512.35

629.735.33.015.4:533.6.013.422

ПРИМЕНЕНИЕ НОВОГО МАТЕМАТИЧЕСКОГО АППАРАТА «ОДНОМЕРНЫЕ СПЕКТРАЛЬНЫЕ ПОРТРЕТЫ МАТРИЦ» К РЕШЕНИЮ ПРОБЛЕМЫ АЭРОУПРУГИХ КОЛЕБАНИЙ

РЕШЕТОК ЛОПАСТЕЙ

В. Г. БУНЬКОВ, С. К. ГОДУНОВ, В. Б. КУРЗИН, М. САДКЕЙН

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

Ключевые слова: качество дихотомии, спектральный портрет, несимметричная матрица.

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

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

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

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

Постулат. Могут вычисляться только такие числовые функции / (А) от N X N или N XМ матриц, для которых справедливо неравенство:

||/(А)-/ (в )|МА - в.

где

известная функция, не зависящая от N, М (порядков матриц);

В|| — нормы матриц.

Пример допустимых функций (а ,■ (А) — сингулярные числа матрицы А):

А = б^Р; дд = 1Ы; РР = Ім; М < М,

0 ^

0

Э =

ам 0

0 ам -1

0 0

0 0

■>м

0 0 0

N.

м

Эта допустимость вытекает из неравенства |а j (А + В) - о j (А) | < а max (В ) = || В11.

Пример недопустимых функций: X}- = Х}- (А) — собственные числа квадратной NXN мат-

рицы А.

Почему это так, можно понять из следующего примера. Взяв целочисленную 7 X 7 матрицу С:

С =

' 289 2044 336 128 80 32 16

1152 30 1312 512 288 128 32

-29 -1980 756 384 1008 224 48

512 128 640 0 640 512 128

1053 2136 -604 -384 -856 800 108

-287 4 1712 -128 1968 -30 2032

ч-2176 -187 -1465 -512 -441 -1152 9 8 7

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

для матриц С и СТ, обязанных иметь совпадающие спектры, различные результаты:

^ (С)= 6.5824 ^ (СТ ) = 8.0444

X 2 (с ) = 4.0313 + 4.3421/ (С) = 4.0313 -4.3421/

X 4 (С ) =-1.4668 + 5.3883/ Х5 (С) = -1.4668 -5.3883/ Х6 (С) = -5.8557 + 2.3387/ Х7 (С) = -5.8557 - 2.3387/

X2 (СТ ) = 4.9557 + 5.6644/ Х3 (СТ ) = 4.9557 -5.6644/

X 4 (СТ ) =-1.8107 + 7.0037/ Х5 (СТ ) = -1.8107 - 7.0037/ Х6 (СТ ) = -7.1672 + 3.0701/ X 7 (СТ) = -7.1672 - 3.0701/.

Матрица С образована по формуле С = Ь 1ЯЬ, где

2028 256 128 64 32 16 > (1 0 0 0 0 0 0'

-2 1024 512 256 128 32 0 1 0 0 0 0 0

0 4 512 1024 256 64 1 0 1 0 0 0 0

0 0 0 512 512 128 , ь = 0 0 0 1 0 0 0

0 0 0 -4 1024 156 0 0 1 0 1 0 0

0 0 0 0 2 2048 1 0 0 0 0 1 0

0 0 0 0 0 -1 V V 0 1 1 0 1 0 1 V

0 0

Я = 0 0 0 0

Откуда видно, что ее собственные значения на самом деле таковы:

X! = 0,

X2 =-1, Xз =+1, X4 =-2, X5 =+2, X6 =-4, X7 =+4.

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

Вот еще один простой пример. Пусть А + юВ матрица 25 X 25, где

О ^

А=

(-1 10 -1 10

О

10

-1

— двухдиагональная, а В имеет лишь один отличный от нуля и рав-

ный единице элемент в нижнем левом углу В =

(0 0

0 0 1 0

—25 —22

При ю = 0 весь спектр А + юВ лежит в левой полуплоскости, а при ю = 10 8 ~2.6 10

среди собственных значений А + юВ существует Х = -4.

После приведенных примеров возникает вопрос: как пользоваться теоремой Ляпунова о том, что для устойчивости решений системы — = А • х необходимо и достаточно, чтобы весь

&

спектр А лежал строго в левой полуплоскости?

Дело в том, что приведенная сейчас формулировка на самом деле — вульгаризация утверждения, доказанного Ляпуновым. Ляпунов показал, что условие Яе(Xу- )< 0 (при всех Xj) необ-

Т

ходимо и достаточно для существования положительно определенного решения Н = Н > 0 матричного уравнения Ляпунова: НА + А*Н + С = 0 при любой положительно определенной С = С* > 0. Построенная с помощью Н, функция Ляпунова (Нх, х) убывает с ростом ^ на решениях системы х = Ах, и поэтому

(' ЖІ Iя

:(°)|.

Каков же коэффициент . Н

(отношение максимальной и минимальной осей эл-

липсоида Ляпунова, в котором лежит траектория х(^))? Этот коэффициент зависит не только

1

от матрицы А, но и от выбора правой части C в матричном уравнении Ляпунова. Таким образом,

*

анализ на устойчивость должен включать не только разрешимость уравнения НА + А Н + C = 0, но и указывать какую-либо конкретную матрицу С, приводящую к приемлемому значению

Я

H

-1

Указанные соображения [1] привели к предложению находить Н из уравнения НА + А Н + 2||А||/ = 0 [1], а характеристикой качества устойчивости считать величину к = ||Н||, при которой справедлива оценка:

-Л

||х(t)||<>/ке к ||х(0)|| = л/ке т ||х(0)||, т = К|А||.

Величина к = к(А) является решением экстремальной задачи:

к(A) = sup J jIх(t)||2 dt jexp(-2^|A||)x(°)||2 dt I.

к

Величина т = — это характерное время для убывания решений.

Неравенство

IBII к-2

справедливое при ^—^.< 10

(А + В )-к(А )< Br’ (A^jjAj

показывает,

что к = к(А) устойчиво по отношению к возмущениям исследуемой матрицы в смысле сформулированного выше постулата.

Выяснилось [2, 3], что Н = Н (А) представляется в виде матричного интеграла:

Я (A ) = U ГГ а* + Ш ] 1 [A -/ю/] 1

d ю,

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

мнимой осью, критерий, оценивающий удаление X}- (А) от этой оси, вне зависимости от того, сколько из точек спектра лежит в правой, а сколько в левой полуплоскости.

График зависимости от а — критерия дихотомии к(А — а1) иллюстрирует расслоение спектра прямыми Яе (X) = а, параллельными мнимой оси. Мы приведем примеры, иллюстрирующие

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

В качестве первого примера, иллюстрирующего использование критерия спектральной дихотомии, рассмотрим простую модель флаттера, предложенную в трудах ЦАГИ [4]. В этой работе крыло-пластина описывается четырьмя степенями свободы. В отсутствие аэродинамических влияний колебания описываются уравнениями:

(37.7 О ^

dx = ^ dy = dt dt

G =

169

899

О

1792

Аэродинамические воздействия моделируются введением в коэффициенты системы новых элементов, зависящих от скорости V потока.

Система превращается в следующую:

■ = ~юх -(+v2 р ),

&у

Б = 0.73 10'

-2

&х

&

&

■ = х,

(1 0^ ( 0 -0.197-10-2 0 0

1 , р= 0.1210-3 0 -0.419 10-2 0.17110-

1 0 0.176 10-3 0 0

V0 1V V 0 -0.154 10-3 0 0

Расслоение спектра этой системы прямыми, параллельными мнимой оси, при различных значениях скорости V иллюстрируется следующими графиками: на рис. 1, а, б — сплошной лиК

нией, пунктирной линией отображен график т = |^ • Удобно совместить отмеченные на них

спектральные зоны на одном рисунке, в котором по оси абсцисс откладывается Яе (X), а по оси ординат — скорость V.

На рис. 1, в серым цветом (внутренняя область, ограниченная штрихпунктиром) заштрихована область значений Яе(X), в которой 1^л/к> 3.95; черная (внешняя область) — при

1§л/К>3.75. Черная средняя область, ограниченная белой линией, — область, в которой К

т =—и>3.75. На рис. 1—6 критерий дихотомии К задается в логарифмическом виде:

ІИІІ

1^л/к = ^ (к)/2, |Х — модуль частоты X, Яе (X) — параметр дихотомии, вещественная часть X. Приведенные графики позволяют оценить суждения о допустимой оценке критической ско-

Яе (Xу)

рости флаттера, получаемой по вычисленному отношению —г—^—

и предлагаемому крите-

рию к(А). Выбор допуска для к(А) должен осуществляться на основе анализа точности моделирования дифференциальными уравнениями и подтверждаться экспериментальными данными.

Отвечающие кластерам собственных значений подпространства (проекторы на них или их базисы) вычисляются одновременно с самими спектральными портретами. Это позволяет указать клеточно-диагональный канонический вид изучаемой матрицы и, сосчитав матрицу подобного преобразования, установить ее число обусловленности. Канонический вид матрицы А (V =411):

(-3.67е + 0 -1.01е + 2 1.17е + 0 -1.32е + 0

в'1 Ав =

0

0

0

0

0

0

0

0

4.30е + 2 3.93е + 2 -2.64е +1 2.68 + 0

-4.76е + 2 -4.33е + 2 2.92е +1 -2.96е + 0

7.85е -1 7.11е -1 -1.79е + 2 1.88е +1

7.75е + 0 7.01е + 0 -1.77е + 3 1.76е + 2

0000 0000

0

0

0

0

0

0

6.73е -1 9.77е +1

-1.21е + 0 -1.67е + 0

в в

1-1

= 267.0132.

Рис. 1. Расслоение спектра прямыми линиями мнимой оси при разных скоростях: а — V=395 м/с; б — V=411 м/с; в — спектральные зоны

В качестве примера для тестовых расчетов спектральных портретов рассмотрим совместные изгибно-крутильные колебания (coupled bending-torsion vibrations) лопаток в решетках турбомашин (blades in the cascades of turbomachinery) в потоке газа.

Система дифференциальных уравнений, описывающих малые колебания такой решетки, имеет вид [5]:

mnhn + Snan + Knhhn = Ln + Fn,

S„h„ + Jnan + Knaan = Mn

(n = 1,2,... N),

где hn, an — обобщенные координаты деформации лопаток при изгибно-крутильных колебаниях; mn, Jn — обобщенные массы и моменты инерции лопаток; Knh, Kna — коэффициенты обобщенных изгибных и крутильных жесткостей; Sn — коэффициенты обобщенной связности изгибных и крутильных колебаний; N — число лопаток в решетке; Fn — силы упругой связности лопастей между собой, вида Fn = c (hn+1 + hn-1 - 2hn); Ln, Mn — обобщенные аэродинамические

силы и моменты, действующие на n-ю лопатку, которые в соответствии с теорией решеток в нестационарном потоке [6] могут быть представлены следующим образом:

N ( h h ' ^

J т/ hr hr I lf \1№ ar

Ln = lr-n,h~T + lr-n,h T + lr-n,aar + lr-n,a ,

, b юЬ ю

r=1 V /

Мп = до^

Г=1

тг—п,к~Г + тг—п,к -Т + тг-п,ааг + тг —п,а ~

Ь юЬ ю

где д — скоростной напор набегающего потока; о, Ь — площадь и хорда лопатки; 1Г-п, тГ-п — аэродинамические коэффициенты влияния лопаток, являющиеся функциями числа Струхаля к = юЬ/У, V — скорость набегающего потока, ю = 1т X, X — корень соответствующего характеристического уравнения этой системы.

Введем обозначения:

юпк ='

К

пк

т„

2 Кпа

Ю = -

^па

п

У п =

А.,

тпЬ'

рп

^,

тпЬ

V =

Ю„

Юпк

е„ =-

до

тпЬю

и отметим, что имеет место оценка £п 1, так как нестационарные аэродинамические силы, дей-

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

+ Упап + ®пл = £п ®21п + — Кп

т„

Упкп + Рп (п +Vю2nкan ) = £пюМп ( = I 2 ... N)

где

N

(

N

г=1 (

\

/' +/" ^г- + /' а +/" ааг.

г-п,к 1 ~ г-п,к 1 ~ г—п,аг г —п,а

Ь юЬ ю

V

;

Г=1

т

г —п,к

- + т

г—п,к~Т + тг —п,ааг + тг —п,а---

юЬ ю

Отметим, что матрица, соответствующая левой части системы, является гамильтоновой, а матрицу, соответствующую правой части, можно рассматривать как ее возмущающую составляющую.

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

юпк, Уп, Рп, V £п, к

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

На рис. 2—6 представлены одномерные спектральные портреты матриц системы, описывающей колебания решетки тонких лопастей при с = 0 (густота решетки т = 1.5; угол выноса в = 30°; прогиб средней линии профиля, отнесенный к его хорде /= 0.025; число Струхаля к =0.5; число профилей в периоде N = 10), обтекаемой идеальной несжимаемой жидкостью. Критерий дихотомии спектра матриц концентрическими окружностями с центром в начале координат комплексной плоскости собственных значений в зависимости от их радиуса на этих рисунках обозначается символом к, а качество дихотомии прямыми, параллельными мнимой оси, от координаты а их пересечения с действительной осью — символом Ха. При этом значения Я и а, при которых к и Ха практически бесконечны, определяют модули и вещественные части собственных значений, а точнее — интервалы для них, за которые можно ручаться на основе приведенных расчетов.

1.5

1.4

1.3

Ні.2 1.1 1.0

0.9

< 3.5

- — 3.0 - 1у ,

\2.5

V мшп

\ 2.0 1.5 - у1 V V щ

/'"Г і і 1.0 .. ,1 . .....

0.5 1.0 1.5 2.0 2.5 -10 -5 0 1

18(к)/2 1ЦХ)-103

а) б)

Рис. 2. Критерий дихотомии к (а) и Ха (б) для решетки лопастей, имеющих одинаковые инерционные и упругие характеристики: у = -0.3, v = 2, р = 1, £ = 0.01

Рис. 3. Параметры дихотомии к (а ) и Ха (б) для решетки лопастей с характеристиками: у = 0, V = 1, р = 1, £ = 0.01

На рис. 2 приведены соответствующие зависимости критерия дихотомии к (рис. 2, а) и Ха

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

у = -0.3, v = 2, р = 1, £ = 0.01.

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

На рис. 3 приведены зависимости качества дихотомии к и Х а для рассматриваемой решетки, характеристики лопастей которой

у = 0, V = 1, р = 1, £ = 0.01.

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

На рис. 4 представлены зависимости качества дихотомии почти для той же решетки, что и на рис. 3, но с учетом упругой связности лопастей между собой. Зависимость качества дихотомии к соответствующей матрицы при с = 0.2юи представлена на рис. 4, а, из которой следует, что данная матрица состоит из жордановых клеток размерности 2 и 4. Согласно имеющимся пред-

Рис. 4. Качество дихотомии для почти такой же решетки, как и на рис. 3, но с добавлением упругих связей в каждой паре лопастей

Рис. 5. Качество дихотомии решетки на рис. 2 (а, б) и изменение его при переходе от матриц большего порядка к меньшему (в, г)

ставлениям такая матрица должна быть более чувствительна к возмущениям, о чем и свидетельствует ее качество дихотомии Ха (рис. 4, б).

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

В качестве примера рассмотрим матрицу, зависимости дихотомий которой представлены на рис. 2. Из зависимости ее качества дихотомии к, с учетом действия упругой связности ¥п (рис. 5, а), следует, что в рассматриваемом случае указанный выше критерий выполняется. Сравнение зависимости качества дихотомии Ха для полной матрицы (рис. 5, б) с аналогичными зависимостями соответствующих подматриц (рис. 5, в, г) подтверждает данное положение.

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

Рп =(-1)П С (п+1 + кп-1 - 2кп ).

Рис. 6. Влияние малых возмущений гамильтоновых составляющих матрицы на примере из рис. 3

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

Описанные в статье приемы исследования устойчивости осуществляются элементарными итерационными алгоритмами, предложенными и описанными в [1, 8, 9]. Алгоритмы решают матричные уравнения Ляпунова и их обобщения на случай дихотомии спектра. Эти обобщения, как нам стало известно, появились в книге [10]. К сожалению, ее содержание не было нами понято, что замедлило обоснование, которому мы уделяли меньше внимания, чем построению расчетных схем.

Авторы благодарят А. А. Сайтгалина за проведение иллюстрационных расчетов.

Работа выполнена при финансовой поддержке заказного интеграционного проекта СО РАН № 5.

ЛИТЕРАТУРА

1. Булгаков А. Я. Эффективно вычисляемый параметр качества устойчивости системы линейных дифференциальных уравнений с постоянными коэффициентами // Сиб. Мат. журнал. 1980. Т. 21, № 3.

2. Булгаков А. Я., Годунов С. К. Круговая дихотомия матричного спектра // Сиб. Мат. журнал. 1988. Т. 29, № 5.

3. Годунов С. К. Лекции по современным аспектам линейной алгебры. — Новосибирск: Научная книга, 2002.

4. Буньков В. Г., Мосунов В. А. Использование интеграла действия по Ляпунову для оценки устойчивости линейной системы // Ученые записки ЦАГИ. 1988. Т. 19, № 2.

5. Bendiksen O., Friedmann P. Coupled bending-torsion flutter in cascade // AIAA. 1980. V. 18, № 2.

6. Горелов Д. Н., Курзин В. Б., Сарен В. Э. Аэродинамика решеток в нестационарном потоке. — Новосибирск: Наука, 1971.

7. Горелов Д. Н., Курзин В. Б., Сарен В. Э. Атлас нестационарных аэродинамических решеток профилей. — Новосибирск: Наука, 1974.

8. Малышев А. Н. Гарантированная точность в спектральных задачах линейной алгебры // Труды Института Математики (АН СССР, Сиб. отд.). 1990. Т. 17.

9. Godunov S. K. and Sadkane M. Some new algorithms for the spectral dichotomy methods of linear algebra // Appl. 2003.

10. Далецкий Ю. А., Крейн М. Г. Устойчивость решений дифференциальных уравнений в банаховом пространстве. — М.: Наука, 1970.

Рукопись поступила 29/V 2009 г.