Аналог метода Хилла в задачах устойчивости периодических решений нелинейных дискретно-континуальных систем Текст научной статьи по специальности «Математика»

ПРОБЛЕМЫ ЕСТЕСТВЕННЫХ НАУК

УДК 531.36

Д.К. Андрейченко

АНАЛОГ МЕТОДА ХИЛЛА В ЗАДАЧАХ УСТОЙЧИВОСТИ ПЕРИОДИЧЕСКИХ РЕШЕНИЙ НЕЛИНЕЙНЫХ ДИСКРЕТНО-КОНТИНУАЛЬНЫХ СИСТЕМ

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

D.K. Andreichenko

THE HILL METHOD ANALOG IN THE STABILITY PROBLEMS OF PERIODIC SOLUTIONS OF NONLINEAR DISCRETE-CONTINUAL SYSTEMS

The analog of the stability Hill method research of periodic solutions with reference to discrete - continual mechanical systems is offered here, driving of which discrete elements are described by the nonlinear ordinary differential equations with lagging argument.

Введение

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

дифференциальных уравнений, в соответствии с теорией Флоке, сводится к вычислению матрицы монодромии и исследованию ее собственных значений [1, с.187].

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

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

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

1. Уравнения движения дискретно-континуальной системы

Уравнения движения показанной на рис. 1 ДКС (совершающего плоское движение абсолютно жесткого спутника с моментом инерции Jc и массой тс под действием возмущающего момента Ь, с жестко закрепленным на расстоянии а от центра массы спутника прямолинейным однородным упругим стержнем, несущим жестко закрепленное на противоположном конце абсолютно жесткое тело с массой т1 и моментом инерции J1) в безразмерной форме имеют вид

Jcа = Ь - /(н>(г -т)) + Ь0 -аЫ0, м>(г) = Ь1<х(г) + Ъ2а(г) + Ъ31 а(^)^Е,, (1.1)

тсУс = N0, Л(« + &&1) = К т1[Ус + - (1+ а)а] = N , (12)

У(z,г) + ус - (а + *)а = -(У""(^ г) + Yy""(z,t)), ( У = д/д* , (1.3)

* = 0 : У (0, г) = У'(0, г) = 0; г = 1: у (1, г) = у (г), у'(1, г) = -а,; (г) , (1.4)

Ьо = - у "(0, г) - ТУ "(0, г), N0 = - у "(0, г) - уу т(0, г)

Ь = у 41, г) + ту41, г), N = у "(1, г) + Уу41, г) ( ' )

г < 0: а(г) = сх(г) = ус(г) = .Ус(г) = у(г, г) = .У(* г) = 0 . (16)

Здесь У(^) - некоторая функция своего аргумента, описывающая нелинейность типа насыщения, Ъ1, Ъ2, Ъ3 - коэффициенты обратных связей ( ) ' = д( )/ дг, ( ) = д( )/ дг. Примем далее

/ ^ = Ш(^0 .

Линеаризация в окрестности состояния равновесия приводит уравнение (1.1) к виду (уравнения (1.2)-(1.6) при этом своего вида не изменяют)

Г1 -т

Jcd = Ь - Ъ1сх(г -т) - Ъ2а(г -т) - Ъ31 а(^)^Е, + Ь0 - а#0, ] = 1,2,3 . (1.1. а)

•< 0

Области устойчивости линейной ДКС (1.1.а), (1.2)-(1.6) в пространстве параметров обратных связей Ъ1, Ъ2, Ъ3, приведенные в [5], подробно исследовались на основе предложенных в [4] теорем об устойчивом квазимногочлене и об устойчивых, неустойчивых и асимптотически устойчивых квазирациональных дробях, а также предложенного в [5] варианта метода .О-разбиений. Следует отметить, что наличие малого запаздывания т резко ограничивало области устойчивости за счет дестабилизации управляемой деформируемой конструкции по высшим формам колебаний. В случае, когда параметры обратных связей выходят за пределы областей устойчивости, исходные уравнения (1. 1)-(1.6) при Ь(г)=0 (с точностью до несущественного выбора начального момента времени) допускают 2п/ш-периодические решения

а = а(0)(гX а = ус = уС°)(гX у = у10)(гX у(г,г) = у(0)(г,г) , (1.7)

которые описывают возможные периодические автоколебательные режимы с частотой ш, исследованные численно в [2] на основе метода гармонического баланса в высших приближениях. Во всех случаях период возможных автоколебаний (даже по высшим формам) превосходил характерное время запаздывания:

2п / ш > т . (1.8)

При численном моделировании развития начальных возмущений [2] выполнялась дискретизация уравнения (1.1) на основе проекционного метода Бубнова-Галеркина с использованием в качестве координатных базисных функций ортогональных полиномов

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

С точностью до выбора несущественного начального момента времени і=0, линейные уравнения движения ДКС, возмущенной в малой окрестности возможного предельного цикла, имеют вид (уравнения (1.2)-(1.5) своего вида не изменяют)

Jcа = -/(0)(і -т)^(і-т) + Ь0 -аЫ0, н,(і) = ^а(і) + Ь2а(і) + Ь3у(і), і>(і) = а(і)

Аналогично [1, с. 298], поскольку исходные уравнения (1.1), (1.2)-(1.5) не зависят от времени, с точностью до несущественного выбора начального момента времени при Ь(г) = 0 уравнения возмущенного движения (2.1), (1.2)-(1.5) по форме записи аналогичны уравнениям относительно набора величин

полученным дифференцированием по времени (1.1), (1.2)-(1.5). Следовательно, 2п/ш-периодические функции

заведомо будут решением уравнения возмущенного движения.

Модельные уравнения (2.1), (1.2)-(1.5) задают некоторый линейный оператор

(строго говоря, задание «предыстории» на отрезке [-т,0] требуется лишь для функций

Далее, поскольку правые части линейных обыкновенных дифференциальных уравнений (2.1), (1.2) являются достаточно гладкими функциями времени, а уравнение (1.3) учитывает малое, но ненулевое вязкое трение в материале конструкции согласно модели Фойгта, относительно негладким «предысториям» У |[-т 0] при г > т будут соответствовать

более гладкие решения У |[1-т г] (рассматриваемые как функции пространственной

переменной и времени). Следовательно, оператор Ц при г > т будет вполне непрерывен в

смысле [6].

Так как коэффициент линейного уравнения (2.1) при слагаемом с запаздывающим аргументом является 2п / ш -периодической функцией времени, при выполнении условия

2. Исследование устойчивости предельных циклов

І(0)(і) = І,(Ь1а(0)(і) + Ь2а(0)(і) + Ь3|а(0)(і)йі), ))г) = 1/еЬ2(г)

(2.1)

Обозначим

У(і) = {у(і), а(і), (х(і), а(і), (х!(іX Ус(і X У с(іX У1(і), У(і X У(г і X У(г і)} .

У(і) = {у(і) = а(іX а(і), с&(і), (х 1(іX а (і), Ус(і), У(іX У(і), У(і), У(г іX У(& і)} ,

(2.3)

у(г) = , а(г), а(г) , для остальных функций требуется задание начальных условий)

(1.7)

(2.4)

и соответственно

Ц+„Г = І, IТ = ІТЦ, Т = 2п / ш, п = 1,2,3,...

(2.5)

Следует отметить, что условие (1.8) не является достаточно жестким ограничением; так, если оно не выполняется, следует выбрать период Т>т, кратный периоду автоколебаний, и соответственно уменьшить частоту ш.

Из (2.5) непосредственно следует, что, в соответствии с теорией Флоке, возможные частные решения уравнений (2.1), (1.2)-(1.5) будут представлять собой произведение показательной функции времени ' на 2п / ш -периодические функции времени

*(>)=I Г««/“. «.(')=е1' 11«“^“. ус о)=е" I;_ ус п п 6)

у()=е- I ;=-; , у(2, ')=г'!;;_; Ус (2у-' .

Далее, легко убедиться, что совокупность значений функций (2.6) при ' е [-т,0] представляет собой «собственный элемент» линейного оператора Ь Т, Т = 2п/ш , а

величина екТ = в2 пХ / ш является соответствующим ему собственным значением. Поскольку оператор Ь Т вполне непрерывен, его спектр ограничен, состоит из не более чем счетного множества собственных значений и имеет единственную возможную точку сгущения вкТ = 0 [6, с. 189], чему соответствует Яе X ^ -; . Таким образом, имеется лишь конечное число решений, которым соответствуют показатели X, лежащие в правой комплексной полуплоскости Яе Х>0, и обеспечивающие неустойчивость предельных циклов. Из (2.2) следует наличие среди (2.6) по крайней мере одного чисто периодического решения, которому соответствует показатель Х=0. Следует отметить, что добавление к показателю X величины гш приводит к сдвигу на 1 индекса Фурье-коэффициентов в (2.6), т.е. фактически подобные показатели неразличимы.

Далее полагаем

Ут(>) = 1 Г=-;Г’е“" • (2.7)

Выражения (2.2) позволяют удовлетворить уравнениям движения (1.2)-(1.5) (в силу их линейности), при этом уравнения относительно ас, а(с1}, у(пс), Уп'>, Уп (2) аналогичны уравнениям (2.2)-(2.5) из [2] при замене Х^Х + шш и приводятся к следующему результату

фу1(Х+ /'пш)ап + фу 2 (Х + /'пш)аС1) +Фv3(X+ гПш)У<п) + (Х+ гСш)2 Фv4(X+ гПш)У(пС) = 0 (2 8)

п = 0, ±1, ±2,..., у = 2,3,4

Подстановка (2.6), (2.7) в (2.1) и перемножение рядов Фурье приводит к следующим линейным уравнениям (вообще говоря, бесконечным)

Ф11(Х + тп)ап +ф12(Х + 1шп)а(1) +ф13(Х + гшп )уП'} + (Х + гшп)2 ф14(Х + тп)у(п0) =

= -(Х + вдс)е-’'Х™п)I У)[*,(Х + /ш/) + Ь2 + Ъs/(X + to/)], п = 0,±1,±2,...

(2.9)

‘ ] *п -

Здесь

ф11(Х) = (^с + а^21 - М-11)Х3 + У0(0)(Ъ1Х2 + Ъ2Х + Ъ3)е ТХ ,

а остальные величины полностью аналогичны [2]

ф14(Х) = Х(«^24 - М"14) , ф12(Х) = (1 + УХ)Х(а^22 - ^12) , ф13(Х) = (1 + УХЖа^23 - ^13) , ф21(Х) = -^21Х , ф22 (Х) = -(1 + УХ)^22 , ф23 (Х) = -(1 + УХ)^23 , ф24(Х) = - М*24 ,

ф31(Х) = (^ 1 - ^31)Х , ф32(Х) = ^1Х - (1 + УХ)Мз2 , ф33(Х) = -(1 + УХ)Мз3 , ф34(Х) = -М'34,

ф41 = -[т1 (1 + а) + д41 ]Х2, ф42 (X) = -(1 + уХ)д42, ф43 (X) = т1Х2 - (1 + уХ)д43, ф44 = т1Х2 - |д44, = ц„ [‘(X)], V,) = 1,2,3,4, ¿(X) = е»'4Л(1 + уЯ)-1/4,

(2.10)

^(/(к) = к- (ch к - cos к - sh к sin к)/а(к), ^(дк) = к- [к (ch к - cos к) - sin к (ch к -1) +

+ sh к (cos к -1)]/а(к), д12(к) = к (sh к - sin к)/а(к), ц13(к) = к 2(ch к - cos к)/а(к) ,

Д14(к) = -ДиЧ к), ^21)(к) = к ~'[sin к (ch к -1) + sh к (cos к -1)]/ а(к),

^21|)(к) = к ~2 [ch к - cos к - к (sh к + sin к) + sh к sin к ]/ а(к), ц22 (к) = -ц13 (к) ,

д 23(к) = -к3(sh к + sin к) / а(к), д 24(к) = -ц^Чк),

д311) (к) = -ДиЧк), д32)(к) = к~3[(1 - ch к)sin к + (cos к + к sin к -1) sh к] / а(к) ,

ц32 (к) = к (ch к sin к - sh к cos к) / а(к), ц33 (к) = к2 sh к sin к / а(к) , ц34 (к) = (к),

= к-2 (cos к - ch к + к ch к sin к + (к cos к - sin к) sh к) / а(к),

|д4?(к) = Д(211)(к), |Д42(к) = |Д33(к), |Д43(к) = к3(ch к sin к + sh к cos к)/ а(к), |Д44(к) = -ц^Чк), а(к) = ch к cos к -1.

Функции ц . (к), V, j = 1,2,3,4 аналитичны в окрестности точки к=0 и при конечных

Y>0 не имеют особенностей в правой половине и на мнимой оси комплексной плоскости (X). Легко проверить

lim X vV u . [k(X)] = b . = const, ReX > 0, v, j = 1,2,3,4

(2.11)

b11 = ß11 =-1/2, b12 = b13 = ^ ß12 =ß13 =0, b14 =-У-1/2, ß14 = 1/2, b21 = (1 + iWlA/Л.

ß21 =-1/4, b22 = b23 = 0, ß22 =ß23 = 0, b24 = л/2у-3/4, ß24 = 3/4, b31 = -/^2у3/4, ß31 =-3/4,

Ьз2 = ->/2y-1/4, ß32 = 1/4, b33 = b34 = y-1/2, ß33 =ß34 = 1/2, b^ = (a + 1)V2y1/4, ß41 =-1/4, b42 =y-1/2, ß42 = 1/2, b43 = b44 = >/2

-3/4, ß43 =ß44 = 3/4 .

Uv. [k (X)] = Uv. [k (X)] .

(2.12)

3. Аналог метода Хилла

Вводя обозначение = ау /(Х + /шу), удобно представить бесконечную

однородную систему линейных уравнений (2.8), (2.9) в форме

1>,+Ф(Х + /шу)е-«Х+”>^/V-0’[Р,(Х + /ш/)2 + Р,(Х + /ш/) + в3]]0, у = 0,±1,±2,... (3.1)

Здесь

Ф(Х) = б(Х)/£>(Х), £>(Х) = ёй{фч.(Х)}, V,] = 1,2,3,4 ; 0(Х) = ёе1{фч-(Х)}, V,] = 2,3,4 (3.2)

Ф(Х), ^(Х), б(Х)- соответственно, аналоги передаточной функции,

характеристического и возмущающего квазимногочленов ДКС [2], [4]. В соответствии с (2.10)-(2.12), функции О(Х) и б(Х) аналитичны при ReХ > 0, и

ImX

©

D(X) = D(X), Q(X) = О (X),

lim X-7D(X) = Const, Re X > 0

(3.3)

Рассуждения, аналогичные п. 2, позволяют

установить, что характеристический квазимногочлен

ReX

Рис. 2

D(X) имеет не более чем конечное число нулей при Re X>0. Заметим, что аналитичность функций D(X) и Q(X), а также конечность числа нулей функции D(X) в правой полуплоскости позволяют достаточно быстро на основании принципа аргумента провести проверку наличия у них совпадающих нулей в правой комплексной полуплоскости. Для этого достаточно сначала проверить наличие нулей D(X) и Q(X) в некоторых вертикальных полубесконечных полосах, а затем, если таковые будут выявлены, проверить наличие нулей D(X) и Q(X) в прямоугольниках некоторой конечной высоты в пределах полуполосы (рис. 2).

Из приведенного в [4] доказательства теоремы об устойчивом квазимногочлене следует, что, в соответствии с (3.3), число нулей (с учетом их кратности) характеристического квазимногочлена D(X) в правой полуплоскости Re X>0 определяется величиной

P = — (7-2 А argD(/о)) . (3.4)

2я V o<a<w /

Будем далее предполагать, что квазимногочлены D(X) и Q(X) не имеют совпадающих нулей при Re X>0; в таком случае величина (3.4) определяет суммарную кратность полюсов передаточной функции Ф(Х) в правой полуплоскости.

Определитель бесконечной системы линейных уравнений (3.1), являющийся аналогом определителя Хилла, вычисляется следующим образом:

D (X) = lim D„ (X), D„ (X) = det{D (X)}, - n <v< n, - n < j < n , (3.5)

n^w

Dvv(X) = 1, D,.(X) = /V:-jV,(“"v^(X + rnv)[b,(X + /fflj)2 + ¿2(X + rnj) + ь,], V] . (3.6)

Используя неравенство Адамара [7, с. 35], можно показать, что в некоторой

конечной области правой комплексной полуплоскости, последовательность (3.5) сходится равномерно, и, следовательно, функция D(X) обладает свойствами, типичными для определителя Хилла:

D(X) = D(X) , (3.7)

D (X + /и) = D (X) , (3.8)

D (X) ^ 1 при Re X^-ro . (3.9)

Сверх того, из равномерной сходимости последовательности (3.5) и свойства периодичности (3.8) следует, что функция D(X) будет аналитична при ReX> 0, за исключением счетного числа полюсов, получающихся сдвигом полюсов передаточной функции Ф^) (т.е. нулей характеристического квазимногочлена D(X)) на величину /пи , n = 0, ±1, ±2,...

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

D (X) = 1 + (D1(X) -1) + (D2(X) - D,(X)) + (D3(X) - D2(X)) +...

значительно улучшается преобразованием в непрерывную дробь [8].

Можно показать, что, если в правой части бесконечной однородной системы уравнений (3.1) заменить нули некоторой последовательностью, для которой сходится сумма квадратов модулей элементов, то решение можно найти по формулам, аналогичным формулам Крамера, и оно также будет представлять последовательность, для которой будет сходиться сумма квадратов модулей элементов. Следовательно, для того, чтобы бесконечная однородная система линейных уравнений (3.1) допускала нетривиальное решение относительно Фурье-коэффициентов vV, v = 0, ±1, ±2,..., необходимо потребовать выполнения условия, аналогичного уравнению Хилла

Аналогично теореме Андронова-Витта можно полагать, что предельный цикл будет орбитально асимптотически устойчив [1, с. 309], если, с точностью до несущественного сдвига на кратную /ш величину, все возможные показатели X будут иметь отрицательные мнимые части, кроме одного показателя Х=0 кратности 1.

Применительно к обыкновенным дифференциальным уравнениям, когда отсутствует запаздывание аргумента, а передаточная функция Ф(Х) представляет собой рациональную дробь, расположение и кратность полюсов которой известны, имеется возможность выделения особенностей определителя Хилла й(Х) и его дальнейшего упрощения на основе теоремы Лиувилля. Однако в случае ДКС сложность поведения передаточной функции Ф(Х) в левой комплексной полуплоскости исключает подобную возможность.

а б в

Рис. 3

Из свойства периодичности (3.8) следует, что, при исследовании устойчивости предельного цикла достаточно исследовать наличие либо отсутствие корней определителя Хилла й(Х) в некоторой полубесконечной полосе (рис. 3, а) в правой комплексной

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

С другой стороны, из (3.6) непосредственно следует, что, если Х=Х0 и Х=Х1 являются различными полюсами передаточной функции Ф(Х), отличающимися на величину, кратную /ш, то Х=Х0 будет полюсом определителя Хилла, причем его кратность будет равна суммарной кратности двух указанных полюсов передаточной функции. Следовательно, суммарная кратность полюсов определителя Хилла в рассматриваемой полубесконечной полосе априорно известна, совпадает с суммарной кратностью полюсов передаточной функции Ф(Х) в правой полуплоскости и определяется величиной (3.4). Однако в подобном случае число нулей в рассматриваемой полуполосе легко устанавливается на основе принципа аргумента.

Значение Х=0 априорно является корнем определителя Хилла й (X), причем

устойчивому периодическому предельному циклу соответствует корень Х=0 кратности 1. Следовательно, с учетом (3.8), (3.9), суммарная кратность полюсов Р и «опасных» нулей N определителя Хилла в полубесконечной полосе (пути обхода контуров 11 и 12 показаны на рис. 3, а и б соответственно) оценивается выражениями

1 Г ^

Р = — 7п-2Ааг§£(Х) I , (3.11)

2п ^ Хе12 у

N = — ilnP- A argD(X))-1 . (3.12)

2n \ Xe11 !

Полуполоса на рис. 3, a должна выбираться так, чтобы ее верхняя и, следовательно, нижняя полубесконечные границы не проходили через полюса определителя Хилла D (X) , т.е. через нули функций D(X + in&), n = 0,±1,±2,.... Такой выбор расположения верхней границы заведомо возможен, т.к. в правой комплексной полуплоскости в любой параллельной действительной оси полубесконечной полосе ширины ш функция D(X) имеет не более чем конечное число полюсов, и может быть достигнут, например, перемещением верхней границы по вертикали (с некоторым шагом) от точки т /2 до точки гш . В свою очередь, проверка отсутствия нулей функций D(X + i^), n = 0,±1,±2,... в узкой окрестности верхней границы проверяется посредством анализа изменения аргумента указанных функций при обходе контура, показанного на рис. 3, в.

В случае, когда входящая в (2.1) функция f (0)(t) мало отличается от константы, т.е. коэффициенты соответствующих модельных уравнений практически не зависят от времени, для коэффициентов Фурье-разложения (2.7) выполняется условие |fv(0)| << fTI.

v = ±1,±2,±3,... В тех областях правой комплексной полуплоскости, где отсутствуют полюса передаточных функций Ф(Х + гпш), n = 0,±1,±2,..., внедиагональные элементы определителя Хилла (3.5) будут пренебрежимо малы. Следовательно, некоторые из возможных нулей определителя Хилла будут содержаться в достаточно малых окрестностях полюсов передаточных функций Ф(Х + шш), n = 0,±1,±2,..., т.е. в малых окрестностях корней (нулей) характеристического квазимногочлена D(X), сдвинутых на величину, кратную гш .

В качестве примера рассмотрим спутник с безразмерными параметрами a=0,05, mc=3, Jc=0,5, m1=1, J1=0,3, y=0,08, т=0,11, ¿1=0,5, ¿2=16,5, b3=5. В данном случае параметры обратных связей не принадлежат области устойчивости, а приведенные в табл. 1

собственные частоты колебаний позволяют сделать вывод о

дестабилизации 3-й формы колебаний деформируемой

конструкции. Система допускает три предельных цикла с

частотами автоколебаний ш=1,25, 4,39, 5,81, показанные на рис. 4, а штриховой линией, пунктиром и сплошной линией соответственно. Фазовые траектории (показаны на рис. 4, б пунктиром) на плоскости (а,а), соответствующие процессу развития малых начальных возмущений, достаточно быстро подходят к предельному циклу с частотой ш=5,81 (показан сплошной линией). Приращение аргумента характеристического определителя D(X) при обходе контура 12, показанного на рис. 3, б, составляет 7п /2, а его годограф в специальном масштабе

u + iv = D(X )ln(1+ | D(X )|)/| D(X )| приведен на рис. 4, в. Следовательно, согласно (3.11), в правой полуплоскости корни характеристического определителя D(X) отсутствуют. Далее, приращение аргумента определителя Хилла D(X) при обходе контура 1Х, показанного на рис. 3, а, равно 2п, а его годограф показан на рис. 4, г. Согласно (3.12), определитель Хилла не имеет «опасных» корней в полосе, показанной на рис. 3, а, и предельный цикл с частотой ш=5,81 устойчив по отношению к малым возмущениям. Ближайшим к мнимой оси корнем определителя Хилла D (X), кроме X = 0 , оказывается X = -0,0368 ± 1,41i.

Таблица 1

Re X Im X

-0,3030 0

-0,0314 1,281

-0,2979 4,600

0,2734 5,989

а

0

-0,2

-0,4

0=1,25 у- I / — 4 о=4,3( \

1/ 7' о=5 ',81 \ п \ \

\ \ \ V #- / 1 У

\ \ \ V Ч.ч. Г

-0,14 -0,07 0 0,07 а

0,15

0

-0,15

-0,3

-0,05 -0,025 0 0,025 а

б

Рис. 4

г

Два других предельных цикла с частотами автоколебаний ш=1,25 и о=4,39 оказываются неустойчивыми. Согласно (3.12), для предельного цикла с частотой о=1,25 в показанной на рис. 3, а полубесконечной полосе определитель Хилла й (X) кроме корня Х=0 имеет еще 3 корня (Х = 0,0561 и Х = 0,1463±0,6002/). Для предельного цикла с частотой о =4,39 определитель Хилла согласно (3.12) имеет 1 «опасный» корень (Х=0,4024).

В качестве другого примера приведем рассмотренный ранее в [2,5] спутник с параметрами а=0,05, дас=34,75, /с=0,07442, да1=3, У1=0,007, у=0,01, т=0,01, Ь1=0,7, ¿2=200, ¿з=150. В данном случае параметры обратных связей не принадлежат области устойчивости, а приведенные в табл. 2 собственные частоты колебаний позволяют сделать вывод о дестабилизации 4-й формы колебаний деформируемой конструкции. Возникает предположение о том, что развитие малых возмущений может привести к возникновению автоколебаний.

Формально, система допускает три предельных цикла с частотами автоколебаний 0=6,787, 11,46, 20,12, показанные на рис. 5, а штриховой линией, пунктиром и сплошной линией соответственно. В начальные моменты времени 0 < I < 4,5 фазовые траектории (показаны на рис. 5, б пунктиром) на плоскости (а,а), соответствующие процессу развития малых начальных возмущений, достаточно быстро подходят к предельному циклу с частотой о=20,12 (показан сплошной линией). Далее, как показано на рис. 5, в, в моменты времени 4,5<^<180 (чему соответствует примерно 560 предполагаемых периодов автоколебаний), фазовые траектории заполняют некоторую относительно узкую окрестность предельного цикла, и становится неясно, является ли предельный цикл

Таблица 2

Ре X 1т X

-0,9948 0

-0,00076 1,003

-0,8905 13,69

-5,550 34,05

6,861 50,55

устойчивым. Приращение аргумента характеристического определителя D(X) при обходе контура 12, показанного на рис. 3, б, составляет -п/2, а его годограф приведен на рис. 5, г. а

-6

ю=6,7 ✓ / \ V |>7 t Г / ч ч

/ t / ю=20, V 12\ \ I

L ' ш= ч у 1,46^ } 1 / /

N > ч г „ + ' /

а

0,7

0

-0,7

-1,4

ю=20, 12 W

А / s:-XV. ■ / / , ! t w » *,v® * v

ь: '• 1 \ * V \ \ \ \ \ Щ: * / < 'S/S. V* 1*’ *♦* - : f ///§ w

|&g * «* .» ,♦

-0,8 -0,4 0 0,4 а

-0,06 -0,03 0 0,03 а

а

3

0

-3

flOl

3 ‘ I At*

и--'""" «

-40

-20

а

д

Рис. 5

а

в

г

0

е

Следовательно, согласно (3.11), О(к) имеет две пары комплексно сопряженных корней в правой полуплоскости (Х = 0,05583 ± 0,9162/ и Х = 0,1251 ± 10,19/ соответственно). Далее, приращение аргумента определителя Хилла й (X) при обходе контура 11, показанного на рис. 3, а, равно нулю, а его годограф показан на рис. 5, д.

Согласно (3.12), определитель Хилла имеет в полосе, показанной на рис. 3, а, две пары комплексно сопряженных корней (Х = 0,04654 ± 0,9373/ и Х = 0,03037 ± 9,136/ соответственно), которые и обеспечивают постепенное отталкивание фазовых траекторий от неустойчивого предельного цикла. Далее, с возрастанием времени ¿, фазовые траектории постепенно уходят на бесконечность (рис. 5, е). Следует отметить, что корни

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

Два других предельных цикла с частотами автоколебаний ш=6,787 и ш=11,46 также оказываются неустойчивыми.

Заключение

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

Работа выполнена при финансовой поддержке Совета по грантам Президента Российской Федерации для государственной поддержки молодых российских ученых и по государственной поддержке ведущих научных школ (грант № МД-2328.2005-1).

ЛИТЕРАТУРА

1. Демидович Б.П. Лекции по математической теории устойчивости / Б.П. Демидович. М.: Изд-во МГУ, 1998. 480 с.

2. Моделирование влияния запаздывающего аргумента в нелинейной газореактивной системе стабилизации спутника с упругим стержнем и закрепленным на его конце телом / К.П. Андрейченко, Д.К. Андрейченко, В.С. Шорин, С.Г. Наумов / Вестник СГТУ. 2005. № 3. С. 17-27.

3. Уиттекер Э.Т. Курс современного анализа. Ч.2. Трансцендентные функции / Э.Т. Уиттекер, Дж.Н. Ватсон. М.: Физматгиз, 1963. 516 с.

4. Андрейченко Д.К. К теории комбинированных динамических систем / Д.К. Андрейченко, К.П. Андрейченко // Известия РАН. Теория и системы управления. 2000. №

3. С. 54-69.

5. Андрейченко Д. К. К теории стабилизации спутников с упругими стержнями / Д.К. Андрейченко, К.П. Андрейченко // Известия РАН. Теория и системы управления. 2004. № 6. С. 149-162.

6. Люстерник Л.А. Краткий курс функционального анализа / Л.А. Люстерник, В.И. Соболев. М.: Высшая школа, 1982. 271 с.

7. Корн Г. Справочник по математике для научных работников и инженеров / Г. Корн, Т. Корн. М.: Наука, 1978. 832 с.

8. Андрейченко Д.К. Эффективный алгоритм численного обращения интегрального преобразования Лапласа / Д.К. Андрейченко // Журнал вычислительной математики и математической физики. 2000. Т. 40. № 7. С. 1030-1044.

Андрейченко Дмитрий Константинович -

доктор физико-математических наук,

профессор кафедры «Механика деформируемого твердого тела» Саратовского государственного технического университета