CC BY
29
УДК 531.36: 62-50 Русанов Вячеслав Анатольевич,
доктор физико-математических наук, главный научный сотрудник, Институт динамики систем и теории управления СО РАН, тел. (3952) 51-17-81, e-mail: [email protected] Антонова Лариса Васильевна, заместитель директора Института математики и информатики, Бурятский государственный университет (БГУ), тел. (33012) 64-89-65, e-mail: [email protected] Козырев Владимир Александрович, аспирант, Иркутский государственный университет путей сообщения,
тел. (3952) 51-17-81, e-mail: [email protected] Шарпинский Дмитрий Юрьевич, научный сотрудник, Институт динамики систем и теории управления СО РАН,
тел. (3952) 51-17-81, e-mail: [email protected]
СТРУКТУРНО-ПАРАМЕТРИЧЕСКАЯ ИДЕНТИФИКАЦИЯ ДИФФЕРЕНЦИАЛЬНЫХ УРАВНЕНИЙ ДИНАМИКИ УПРУГОГО ЭЛЕМЕНТА
СПУТНИКА-ГИРОСТАТА
V.A. Rusanov, L. V. Antonova, V.A. Kozyrev, D. Y. Sharpinsky
STRUCTURAL AND PARAMETRIC IDENTIFICATION OF THE DIFFERENTIAL EQUATIONS OF THE DYNAMICS OF THE SATELLITE-GYROSTAT ELASTIC ELEMENT
Аннотация. В работе представлен программно-алгоритмический комплекс, имитирующий этапы математического проектирования измерительной аппаратуры и вычислителя, в апостериорном моделировании уравнений многомерной дифференциальной системы, описывающей динамику демпфированных поперечных колебаний защепленной штанги-балки спутника-гиростата. Функциональной основой данных этапов является формализованная (в рамках метода Фурье) аналитическая модель изгибных колебаний балки Эйлера - Бернули, представленной системой дифференциальных уравнений счетного пучка стоячих волн. Предложен метод идентификации, дающий возможность диагностирования элементов крупногабаритных стержневых конструкций.
Ключевые слова: спутник-гиростат, уравнения движения в пространстве, дифференциальная реализация.
Abstract. This paper presents a software-algorithmic complex to simulate the stages of mathematical design of the measuring apparatus and a calculator in aposterior simulation of multidimensional differential equations system describing the dynamics
of damped transverse vibrations of satellite-gyrostat beam. Functional basis of these stages is a formalized (in the Fourier method) analytical model of flexural vibrations of Euler - Bernoulli beam, provided with the system of standing waves countable beam differential equations. Method of identification, enabling diagnosis of large elements of the bar structures, is proposed.
Keywords: satellite-gyrostat, the equations of motion in space, a differential implementation.
Введение
Голономную механическую систему (голо-номная механическая система - это система, в которой все механические связи являются геометрическими и/или налагающими ограничения только на положения (перемещения) точек и тел системы, но не на их скорости, как это имеет место в него-лономных системах [2, с. 13]), содержащую в качестве своих частей как подсистемы с конечным числом степеней свободы, так и звенья с распределенными параметрами, для краткости, как правило [1], называют сложной механической системой (СМС) (такая «негомогенная» система зани-
Современные технологии. Механика и машиностроение
ш
мает промежуточное место в последовательности возможных механических моделей между одиночным абсолютно твердым телом и сплошной средой; в зарубежной литературе подобную систему уравнений иногда называют гибридной [3]). В такой постановке последняя есть механическая система со счетным числом степеней свободы. Соответственно, ее движение описывается совокупностью счетного числа обыкновенных дифференциальных уравнений, в общем случае нелинейных и с переменными (нестационарными) коэффициентами. Спутник-гиростат (СГ), содержащий в своем составе упругий элемент типа гибкого стержня (твердый гиростат + гибкий стержень) - важный пример сложной голономной механической системы:
Дифференциальные уравнения СГ-динамики структурно и параметрически определяются геометрией СГ-системы, распределением масс в СГ-системе, а также природой внешних сил и сил, действующих в местах соединения тел СГ-системы. При пространственно-угловом СГ-движении эти взаимосвязи эволюционируют, причем по законам, не поддающимся априорному описанию, что выдвигает задачу их апостериорной оценки, а если говорить шире - структурно-параметрической идентификации нелинейных уравнений СГ-динамики (хотя основные результаты этой работы касаются именно модели спутника-гиростата, будет ясно, что многие рассмотрения обратных задач динамики можно перенести на случи более общие, чем СГ-динамика).
С появлением космических СМС необходимость в учете влияния нежесткости их конструкций становится все более насущной; наличие в геометрической конфигурации СМС деформируемых звеньев может существенно изменить динамические характеристики, в частности определяющие устойчивость орбитальной ориентации системы (ее стационарных положений). Эту пробле-
му решает качественная теория систем. Важный аспект теории нелинейных динамических систем -исследование и классификация точек равновесия в зависимости от структуры и параметров модели исследуемой системы, в том числе при ее адаптации к условиям реальной технической эксплуатации.
Как показывает практика, проведение точного количественного анализа нелинейной динамики СМС возможно только в редких случаях: среди нелинейных дифференциальных уравнений, описывающих динамику сложных систем, лишь для немногих удается получить точное аналитическое решение. Однако часто на практике не только важно знать точное численное значение параметров системы, задающее (определяющее) ее состояние в заданные моменты времени, но и описание поведения системы в целом, т. е. на качественном уровне. Иными словами, при проектировании и текущем поддержании функционирования СМС важно знать, существуют ли у нее состояния равновесия, какие из них устойчивы, а какие - нет, возможны ли в системе устойчивые колебательные режимы и т. п. Необходимо знать также, при каких условиях в системе возможны те или иные характерные поведения (колебания, равновесия и т. п.) и когда устойчивые состояния сменяются неустойчивыми.
Сказанное выше делает актуальным структурно-параметрическое моделирование динамики сложной системы [4], как при ее стендовых испытаниях (или испытаниях отдельных подсистем), так и в режиме ее длительного функционирования. Поэтому ниже основная целевая установка работы - описание алгоритмов комплекса «РЕ-ДИМ-МОДФОКБ» по дифференциальной реализации экспериментальных данных для гибкого аксиального элемента спутника-гиростата в виде нежесткой штанги.
1. Постановка задачи апостериорного моделирования динамики голономной системы с большим числом степеней свободы
Рассмотрим СМС, состоящую из конечного числа твердых, упругих, а также материальных точек (т. е. с ненулевой массой). Одно из твердых тел СМС примем за основное (в системе связанных тел обычно выделяют одно основное тело (как правило, наиболее массивное); А.И. Лурье назвал [6] такое тело носителем системы, а все остальные тела, с ним связанные, - носимыми) - p0. С телом p0 свяжем одно упругое тело плотностью р (случай нескольких упругих тел легко получается из рассматриваемого ниже) и некоторых других твердых тел и материальных точек; полости с жидкостью и т. п. не рассматриваем.
С носителем p0 жестко свяжем связанную прямоугольную декартову систему координат ОХ1Х2Х3 с началом в центре масс O тела p0 и осями, совмещенными с его главными центральными осями инерции. Абсолютное движение системы будем рассматривать по отношению к инерциаль-ной системе координат О начало которой,
например, помещено в центре планеты, а оси направлены так, чтобы в начальный момент времени они были параллельны осям орбитальной системы координат; начало орбитальной системы помещается в характерную точку СМС, например O, одна ось направляется по радиусу-вектору, другая - по нормали к плоскости орбиты, третья - по трансверсали [5].
Положение какой-либо точки py системы определим ее радиусом-вектором
Гу'=Го'+Гу
относительно точки О ', где г0 ' - радиус-вектор точки O относительно О' , гу - радиус-вектор точки Py относительно O (в координатах осей Oxt, i=l,...,3).
В такой постановке абсолютная скорость точки py равна:
vy=v0+<xry+dry/dt, где x - операция векторного произведения, v0 и ю - векторы скорости точки O и мгновенной угловой скорости тела p0, dry/dt - вектор скорости точки py в ее движении по отношении к телу p0.
Кинетическая энергия системы, как известно [6, 7], имеет следующий вид: E=mEv02/2+mEv0T(raxrc)+ra т©0ю /2+v0TEymydry/dt+ +ю Eyryxmydry/dt+Eymy (dry/dt)2/2, где mE=Eymy - общая масса СМС, гс - радиус-вектор центра масс СМС относительно точки O, ©0 - тензор инерции системы для точки O; знак Ey означает суммирование по всем точкам СМС с массами my (здесь и ниже части сумм Ey, относя-
щиеся к частицам упругого тела, представляются интегралами, взятыми по области D пространства Oxix2x3, занятой в данный момент времени упругим телом), «т» - операция транспонирования вектора.
Количество движения Q и момент количеств движения G относительно точки O определяются по формулам (где vc - вектор скорости центра масс СМС):
Q=gradv0, E=mEvc, G=gradffl E=rc x mEv0+©0< +Eyryx mydry/dt.
В результате будем иметь замкнутую систему нелинейных дифференциальных уравнений СМС - обыкновенных и в частных производных [6-8]:
dQ/dt+<axQ=F, dG/dt+<axG+v0xQ=M, pd2w/dt2=fdw/dt-Lw+U, (1)
Lw:=-div(p(x)grad w(x,t))+q(x)w(x,t), xeD, p(x)d2 w(x, t)/dt2=fx)dw(x, t)/dt+div(p(x)grad w (x, t))--q(x)w(x,t)+U(x,t), где F и M - главные вектор и момент относительно точки O всех приложенных к системе активных сил Fy, величина w(x,t) означает относительное перемещение точки py при деформации упругого тела (в частности, это относительное перемещение суть любая из координат вектора ry), функции p(x), fx), p(x) и q(x) определяются свойствами среды, где происходит деформация упругого тела, а свободный член U(x,t) выражает интенсивность внешних воздействий Fy.
Распространенный способ дискретизации СМС состоит [8, с. 464-497] в представлении перемещений упругих элементов системы в виде рядов по некоторой полной системе пространственных функций (по собственным функциям эллиптического оператора), умноженных на обобщенные координаты, при этом в указанных рядах оставляют небольшое число первых членов. Этот метод обычно называют «методом нормальных форм колебаний» (менее часто - «методом усечения числа гармоник»). Один из недостатков метода нормальных форм колебаний состоит в том, что усечение указанных выше рядов вносит элемент неопределенности (например, в условиях устойчивости равновесных положений СМС), поскольку затруднительно заранее определить число форм колебаний, необходимое для того, чтобы этот метод дал удовлетворительные результат. В свою очередь, отмеченное положение делает актуальным экспериментальное подтверждение (или опровержение) «рекомендаций метода усечения», что формализует следующая задача апостериорно-
го моделирования СМС (модификация задачи из
[9]).
Выделим класс систем [9], движение кото" nn
рых в пространстве состояний R описывается век-торно-матричным дифференциальным уравнением dx(t)/dt=Ax(t)+By(x(t)), teT:=[t0,t!], (2) x(to)=XoeRn, y(t)=Cx(t);
xeAC(T, Rn) - траектория системы, y(x)eL2(T, ц, Rm) - нелинейная составляющая уравнений состояния системы (или управление в форме нелинейной обратной связи), yeAC(T,Rp) - «выход» системы, AeMnn(R), BeMn,m(R), CeM^R), p<n, Mn_m(R) — пространство всех nxm-матриц с элементами из R.
Задачу дифференциальной реализации [10] будем рассматривать (с учетом второго уравнения в (2)) в рамках следующего структурного ограничения
V(x(t))=9(y(t)), teT, (3)
обусловленного фактом: гипотезу закона y(x) можно экспериментально подтвердить, или опровергнуть исключительно с помощью выходного сигнала y(t).
Постановка задачи структурно-параметрической идентификации: для апостериорного процесса yeAC(T,Rp) и априори заданного (как гипотеза) закона ф: AC(T,Rp)^L2(T,^,Rm) таких, что
dim Span{yeAC(T,R): i=1,...,p}=p,
dim Span^;(y)eL2(T,^,,R): j=1,...,m}=m, (4) Span^;(y): j=1,...,m}, Vi=1,...,p;
col(y1(t),^^yp(t)):=y(t), col^(y(t)), .. ,^m(y(t)))^(y(t)), получить конструктивные процедуры в решении следующих задач:
1) определить необходимые и достаточные условия разрешимости задачи дифференциальной реализации (совместное выполнение позиций а) -в), что по существу суть структурная идентификация моделируемой системы):
а) существует конечномерное фазовое многообразие Rn минимального динамического порядка n (называемое минимальным пространством состояний); ясно, что размерность n минимального пространства состояний опосредованно определяет число нормальных форм колебаний при построении редуцированной модели СМС в рамках метода усечения гармоник,
б) существует некоторое начальное состояние x(t0)=x0eRn,
в) существует дифференциальная система (2), (3) (с некоторыми матрицами A, B, C), эволюционирующая (с начальным условием x0) на фазо-
вом многообразии Я" и имеющая такое же (идентичное) отображение у=Сх;
2) в положении 1) построить прямой алгоритм реализации (параметрическая идентификация системы): вычисление минимального динамического порядка п, вектора х0еЯ", матриц ЛеМпп(Я), БеМпт(Я), СеМрп(Я); в общем случае конструкции х0, Л, Б, С далеко не единственны.
Первое условие в (4) по существу означает, что любая переменная у^ наблюдаемого (экспериментально) на интервале Т векторного выходного СМС-сигнала у не может быть выражена посредством линейной комбинацией других у}- (это приводит к минимальной размерности пространства выходных сигналов); из второго условия (4) следует аналогичное структурное свойство для координат вектор-функции ф(у). Третье условие специфично, в определенном смысле оно констатирует: реализация должна быть с пространством состояний, поскольку из него с необходимостью следует:
уфЮф), У£еМр,т(Я).
Пункт а) позволяет называть такую реализацию минимальной [10, с. 267]; в обычных матричных терминах этот пункт эквивалентен положению [10, с. 68], когда пара матриц (С, Л) является наблюдаемой (или, для метода нормальных форм колебаний, - наблюдаем усеченный ряд гармоник).
2. Модели внешнего упругого элемента спутника-гиростата в форме уравнений Ла-гранжа второго рода
Движение сложной голономной механической системы (как и движение любой механической системы) удовлетворяет основным принципам и теоремам аналитической механики, отраженным в дифференциальной системе (1). Так как любая СГ-система есть частный вид сложной механической системы, целесообразно преобразовать общие СМС-уравнения (1) к более простому (частному) выражению, соответствующему именно дифференциальным уравнениям (таким дифференциальным уравнениям отвечают уравнения Ла-гранжа второго рода, поскольку система уравнений в независимых координатах имеет [2, с. 50] наименьший возможный порядок, что отвечает а) п. 1.) СГ-динамики; поскольку аналитические трудности нелинейной теории упругости довольно значительны, ее приложения к вращающимся СГ-конструкциям, как правило [11], ограничиваются случаями радиальных и аксиальных стержней, а также осесимметричной пластиной, перпендикулярной оси вращения спутника-гиростата (или даже ограничиваются одной из данных схем -рис. 1).
В такой постановке анализ (в рамках метода усечения числа гармоник) форм колебаний упругого элемента (в частности гибкого стержня) на вращающемся СГ-корпусе требует вывода квазилинейных уравнений (2), описывающих малые колебательные отклонения элемента от его постоянного деформированного состояния [12]. При этом структура гипотезы нелинейного закона (3) должна определяться методологической установкой, что основной фактор нелинейности СГ-динамики по существу вызван вращением носителя СГ-системы, т. е. телом p0.
Рассмотрим аналитическую модель уравнений СГ-динамики в форме дифференциальных уравнений Лагранжа второго рода [11, с. 101]: Md2w(t)/dt2+Ddw(t)/dt+Sw(t)=Ku(w(t),dw(t)/dt), (5) w(t0)=w0, dw(t)/dt | t=t0=w'0, y(f)=Cco\(w(f),dw(t)/df), где wе - вектор-функция обобщенных
перемещений (для СГ-конструкции вектор-функция д=^+д2 декомпозируется на вектор-
функцию ^ - перемещения конструкции как «жесткого целого» + упругих перемещений - д2.), иеС <ю(T,Rm) - вектор-функция внешних нагрузок, уеС Ш(Т,ЛР) - измерения датчиков первичной информации, Ми S - симметричные кхк-матрицы масс и жесткостей конструкции, D - кхк-матрица демпфирования, К - кхт-матрица внешних воздействий, С - рх2к-матрица наблюдающего устройства, w0, w '0 - начальные векторы перемещений и скоростей (состояние конструкции в момент времени
З а м е ч а н и е 1. Система (5) записана в предположении о малости вектора перемещений; это накладывает определенные ограничения на
вид вектор-функции f^u(w(t),dw(t)/dt) внешних
нагрузок [13, с. 93]. Поэтому применительно к СГ-динамике дифференциальные уравнения (5) могут использоваться на ограниченном интервале времени при рассмотрении случаев приложения к конструкции уравновешенных [14, с. 275] и слабо неуравновешенных воздействий, когда движение конструкции как твердого целого мало отличается от равномерного поступательного движения [11, с. 102].
Уравнения Лагранжа (5) легко трансформируются к виду Коши (2) (преобразование к модели в пространстве состояний) через линейную связь: х(0:=со\М0, dw(f)/df)еRn, п=2к, (6) у(х(ф:=и^(^, dw(f)/df),
А :=
0,
'к Ек - М -1S - М 1D
В :=
0к М 1К
где 0к и Ек - соответственно нулевая и единичная кхк-матрицы. Ясно, что представления (6) приводят математическую постановку контекста задачи апостериорного моделирования из п. 1 к задаче структурно-параметрической идентификации квазилинейных уравнений СГ-динамики в форме (2).
3. Разработка интерактивного комплекса «РЕДИМ-МОДФОКБ» для моделирования динамики поперечных колебаний гибкого аксиального элемента спутника-гиростата
Выше были определены условия структурно-параметрической идентификации распределенной динамической системы, отвечающей схеме метода Фурье (представление решений распределенной системы в виде сумм стоячих волн), в частности для поперечных деформаций «длинной» аксиальной штанги-балки, жестко прикрепленной к спутнику-гиростату (рис. 1). Для этой цели служит программная среда «РЕДИМ». Программная среда «РЕДИМ» - это живой организм, она постоянно развивается: в функциональное наполнение добавляются новые модули, старые модифицируются; разумеется, каждое вносимое в тот или иной модуль изменение должно автоматически отразиться во всех существующих в данный момент копиях этого модуля - и тех, которые хранятся на внешних носителях, и тех, которые участвуют в сборке отдельных программ. В работе [15] был разработан интерактивный комплекс «МОД-ФОКБ» (моделирование динамики форм-колебаний балки), основанный на методе усечения числа гармоник-деформаций балки с жестко закрепленным концом. Графический интерфейс в программной среде «МОДФОКБ» выполняет функции сервисного характера [22], и основные задачи здесь состоят в следующем.
Прежде всего, необходимо организовать хранение функционального наполнения. Функциональное наполнение представляет собой набор отдельных программ, решающих конкретные задачи. Эти задачи объединены общей направленностью, или, как говорят, предметной областью. Дело в том, что программная среда «РЕДИМ» не является универсальной, она проблемно ориентирована, т. е. предназначена для решения определенного класса задач апостериорного математического моделирования динамических систем. Но хранить не значит только ограничиться записью информации на каких-либо носителях; в архиве должен быть порядок (по первому требованию указанный модуль должен быть включен в вычислительную работу). Поэтому второе назначение «МОДФОКБ» - обеспечить возможность сборки из отдельных модулей полной программы, спо-
собной решать задачу моделирования. Для этого интерфейс, создающий граф решения предметной задачи, должен общаться как с пакетом «РЕДИМ», так и с пользователем: давать заказы, воспроизводить ответную информацию и т. п. Этой цели служит язык заданий, т. е. средство диалога пользователя с пакетом, реализованное на языке Delphi. При этом пользователь разбивает решаемую предметную задачу на части, подбирает для решения каждого куска свой проблемно-ориентированный модуль, и затем «МОДФОКБ» начинает процесс соединения модулей в единый вычислительный комплекс; это упрощенная схема, но она отражает характерные этапы работы. Модульная структура «МОДФОКБ» выполнена с использованием среды MATLAB (пакет «Control System Toolbox») и обеспечивает решение следующих задач:
Модуль-1: формирование «банка» стоячих демпфированных волн-гармоник, описывающих поперечные колебания гибкой штанги с защепленным концом.
В Модуле-1 формирование собственных форм-колебаний wk стоячей волны поперечного колебания штанги можно проводить по усмотрению пользователя, либо в виде сплайн-аппроксимации числовых массивов, задающих данные формы (см. ниже задачу (13)-(16) моделирования балки Тимошенко), либо (см. замечание 2 к моделированию балки Эйлера - Бернули) в аналитическом виде
x\^wk(x)=vk-(cos(0,5k%l'lx) - 1), vkeR, xe[0, l]; обоснование полноты системы {wk}k=1, 2,... см. в [17, с. 393]). В такой постановке k-я стоячая волна yk демпфированного поперечного колебания имеет вид
(x,t)\^yk(x,t)=wk(x)exp(-bkt)sin(akt+^k), (x,t)e[0, l]xT, (7)
(k = 1 - основной тон, k Ф 1 - обертон) с моделью датчика, измеряющего в аксиальной точке x0e(0, /] текущую амплитуду «прогиба» волны yk:
t^Yk(xo,0=Vk-(cos(0,5krc/-1Xo) - 1)exp(-bkt)sin(<Bkt+9k), VkeR, teT.
Функциональные характеристики k-й стоячей поперечной волны
T=[0, т] - интервал моделирования (по умолчанию х=10 сек);
l - длина штанги (/=10 м, x0=0,7/); (ak - собственная частота k-й формы колебаний (0,5Ы-1);
q>k - сдвиг фазы k-го обертона, 0<9k<rc (9k=0,5n);
vk - коэффициент масштабирования амплитуды k-й волны (vk=k~2);
bk - коэффициент затухания k-й гармоники
(bk=0,1k0'3).
В Модуле-1 функционально предусмотрено:
1) имитационное моделирование демпфированного колебания изолированной (отдельной) k-й поперечной стоячей волны (x,t)—yk(x,t) штанги;
2) сплайн-интерполяция вдоль линии профиля штанги собственной формы x—wk(x) для k-й
гармоники осуществляется с шагом Ax=0,02k_1l;
3) на временном интервале T анимация динамического процесса колебания каждой k-й стоячей волны проходит с шагом At=0,01k_1T сек).
Графическая часть Модуля-1 позволяет (рис. 2):
1) строить на отрезке [0, /] форму распределенной амплитуды xi—¥wk(xx) (осуществляется методом кубической сплайновой аппроксимации);
2) строить анимированный на временном интервале T динамический процесс (x,t)—yk(x,t)
РЕДИМ - МОДФОКБ
Рис. 2. Фрагмент экранной формы интерфейса Модуля-1 (вариант k=5)
колебаний к-ой стоячей волны в профиле штанги;
3) строить анимированный на Т сигнал ^ук(х0,0 датчика и его производной
f^dyk(xo,t)/df=-bkУk(xo,t)+ЮkWk(xo)e~шcos(Юkt+фk) в аксиальной точке х0е(0, /].
З а м е ч а н и е 2. Использованное в Моду-ле-1 аналитическое представление собственных изгибных форм wk можно назвать «упрощенная модель балки Эйлера - Бернули», поскольку задача (13)-(16) в постановке г2=д2=0 (см. ниже замечание 6) приводит к дифференциальному уравнению
diwk(x)/dx4 - Хк^к(х) = 0 с частным решением (для защепленной балки Эй-лера-Бернули) вида:
xl^wk(x)=vk•(cos(A,k0'5x) - сЬ(Як0'5х)), хе [0, I], к=1,2,...
Модуль-2: моделирование на базе принципа суперпозиции динамики кортежа поперечных форм-колебаний напряженно-деформированной штанги со свободным торцом (сборка функционального пакета из N различных волн-гармоник).
В Модуле-2 функционально предусмотрено:
1) задав целое N проводить фиксацию {..., .. различных к-номеров волн-гармоник ук из Модуля-1 для сборки функционального пакета ^ ((1-4,6,7Ь);
2) для каждой гармоники из выбранного пакета ^ определять и обозначать графически серии точек для узлов {..., ...}к.узлы и пучностей
{. • ^ .}к-пучности;
3) задавать число р измерительных датчиков для точечной съемки информации об амплитуде моделируемого колебательного движения штанги (р = N);
4) выбирать на штанге точки с аксиальными
координатами {хь ..., хр}с(0, /], 0<хг</, 7=1,., р, расположения датчиков съемки информации (р=^ см. ниже утверждение 1), измеряющих текущий «суммарный прогиб» демпфированных волн-гармоник функционального пакета
({ . • ^ .}к-пучности).
Аналитическая модель функционального пакета ^ N-усеченного ряда числа гармоник свободных колебаний гибкой штанги в Модуле-2 имеет вид
(х,0^(х,0=11,...лУк(х,0, (х,0е[0, Г]хТ, (8)
при этом модель измерений датчиков <^(х,01 х=х7, 7=1,., р, представлена как
^^(хь0=Е1,...лУк(хь0, tеT;
^^(хр,0=!1,...лУк(хр,0, tеT. Графическая часть Модуля-2 позволяет:
1) строить (сплайновая интерполяция линии изгиба штанги х^^^хд)^ в любой момент t проводится на равномерной сетке с шагом Ах=0,02ктах_1/, при этом анимация на временном интервале Т процесса колебания «динамического пакета гармоник» осуществляется с постоянным шагом Аf (А?=0,01 с), анимированное на интервале Т колебание штанги (х,0^^(х,0, представленное
функциональным пакетом ^усеченного ряда числа гармоник (рис. 3);
2) строить анимированную на интервале Т функцию х7е(0, Г] сигнала 7-го (7=1,.,р)
датчика и его производной f^dЪlN(xi,t)/dt (рис. 3);
3) строить анимированную на временном интервале Т функцию линейной свертки
^(Е1,. ,р д£4хг,0+Е1,. ,р g1dl¡N(Xг,t)/dt), дг, gгеR от
Рис. 3. Фрагмент экранной формы интерфейса Модуля-2 (вариант {1,5}м=2)
сигналов датчиков поперечных колебаний и их производных.
Модуль-3: диагностирование частот волн-гармоник поперечных деформаций гибкой штанги посредством апостериорной оценки спектральной плотности мощности сигналов датчиков величины ее прогиба (преобразование Фурье автокорреляционных функций 1^^()=СХ}\'ЮЪ^(хьт)Ъ^(хьТ+1)ёт, 7=1,...,р).
В Модуле-3 функционально предусмотрено:
1) проводить спектральный анализ модели посредством преобразований Фурье автокорреляционных функций 1^^г(1)=ш.Рю^(хьт)^(хьт+1)^т,
¿=1,., р, методами Томсона, Уэлча и периодограммы с помощью МЛТЬЛБ-операторов 1)1111:111, pwelch, periodogram;
2) выделять для сигналов й^^(х7,1), хге (0, /],
¿=1,., р, частоты их максимальной спектральной плотности с целью: оценки числа N активированных в штанге форм-колебаний, оценки значений собственных частот ю к деформаций штанги, оценки коэффициента \к масштабирования амплитуды к-й стоячей волны.
Графическая часть Модуля-3 позволяет:
1) выполнить МЛТЬЛБ-процедурами рийт, pwelch, periodogram оценку (понятие оценки (приближенного значения) весьма важно, т. к. на практике почти исключительно оперируем приближенными величинами; это происходит либо потому, что точные значения исходных данных неизвестны (измерение всегда производится неточно, с большей или меньшей погрешностью, в зависимости от точности измеряющего инструмента), либо потому (наш случай), что действия,
которые следует выполнить над данными, были бы слишком сложны и кропотливы, если искать «точные значения» искомых величин; таким образом, в расчетах пользуются тем или иным приближением, в зависимости от того, какая точность нас устраивает) спектральной плотности | £7(ю) 12=<х1'ю^г(т)е"-'(ЮУт мощности (спектральная плотность сигнала 1^^(хь1) суть комплексная величина 57(ю )=<арю^м(х7,т)е^ш'Ут, поэтому спектральную плотность мощности сигнала определяет преобразование Фурье его автокорреляционной функции) сигналов 11^^(хь1), ¿=1, ..., р [38, с.74],
строя их графики в логарифмической шкале (ремарка: данная операция эффективна при больших частотах юк (см. рис. 4), либо на малых частотах юк при невысоком («0) коэффициенте затухания Ьк);
2) делать оценку частот сигнальных функций [18] для алгоритма параметрической идентификации, реализованного в Модуле-4.
З а м е ч а н и е 3. В идеале программный комплекс должен быть оснащен инструментом, который автоматически решает задачу пункта 2, однако это перспективная цель, а пока пользователь на эмпирической базе, полученной средствами рт:т, pwelch, periodogram, сам должен заниматься данной оценкой.
Модул ь-4: структурно-параметрическая идентификация дифференциальных уравнений поперечных колебаний гибкой штанги с жестко закрепленным концом (дифференциальная реализации сигналов датчиков 1^^(хь1), 7=1,...,р).
В Модуле-4 функционально предусмотрено:
1) определение количества всех активированных в Модуле-2 осцилляторов волн-гармоник
1—
^ Моделирование поперечных колебаний гибкого элемента спутника-гиростата Версия 2.0.0.2
Рис. 4. Фрагмент экранной формы интерфейса Модуля-3 (вариант {1,5}м=2)
пакета (x,t)—^N(x,t) на основе выбора числа и геометрии размещения датчиков измерений t—^N(xbt), i=1,...,p (см. ниже утверждение 1);
2) построение матрицы А дифференциальной модели, описывающей динамику системы - «свободные колебания гибкой штанги (x,t)—^N(x,t) &
сигналы измерительных датчиков t—ílN(xi,t),
i=1,.. .,N» (в силу утверждения 1):
dx(t)/dt=Ax(t), teT, x(t)eR2N, (9)
x(t):=col(x1(t),.,X2N(t)),
fxi(t)=^N(xi,t), I X2(t)=d^N(xi,t)/dt,
í.....................
1 X2N-1 (t)=^N(xN,t), lx2N(t)=d^N(xN,t)/dt
с привлечением алгоритма из [9] в построении линейной однородной системы дифференциальной реализации апостериорных данных
t—(^N(xi,t),dí¡N(xi,t)/dt), i=1,...,N (при полном измерении вектора состояния col(x1(t),.,x2N(t)));
3) обработка спектра оператора системы (9) на базе MATLAB-операторов: esort (сортировка комплексных собственных значений системы), damp (вычисление собственных частот и коэффициентов демпфирования для полюсов системы), pzmap (расположение полюсов и нулей);
4) приведение (оператор canon) матрицы A к канонической форме Фробениуса с вычислением матрицы C наблюдающего устройства y^^x^).
З а м е ч а н и е 4. Алгоритмическое решение задачи 1 Модуля-4 достигается на основе решения двух задач: определения числа N и выбора мест установки датчиков t—^N(xbt), i=1,...,N что формализует утверждение 1.
У т в е р ж д е н и е 1. Для пакета волн-гармоник (x,t)i—^N(x,t) в аксиальном отрезке [0, Г]
существует набор из N точек {xb...,xN}c(0, /], такой, что для измерительных датчиков t—^N(xi,t), teT, i=1,...,N, выполняется отношение:
r^N^!,-), ..., £n(xn,0, d^N(xb•)/dt, ..., dUxN,-)/¿t)*0, (10)
при этом для любой точки ze(0, /] и датчика t—^N(z,t), teT имеет место
шшт
^(^(хь-), ..., Ы%,-), Ш^Ь^М, ..., Ш^^ОШ ^(г,0)=0. (11) В данных условиях пакет ^ имеет дифференциальную реализацию (23).
Д о к а з а т е л ь с т в о. Доказательство проведем индукцией по числу N.
Базис индукции N=1. В силу (7) Е,1(х,^=ук(х^^к(х)еЬк^т(а^+фк), хе[0, I], tеT. Пусть х1е(0, /] - узловая точка волны ук+1, такая, что ^ук(хь0^0. Тогда
^1(х1,0=се-ьк^т(ю ^+фк), d^1(x1,f)/df= =-сbkЮke"bktcos(ю kt+фk), с=Wk(Xl )*0.
Прибегая к теореме 1 [14, с. 213], заключаем, что справедливы условия (10), (11) и, значит, для утверждение 1 осуществимо с {х1}.
Шаг индукции от N-1 к N. В силу (8) ^(х,0=^-1(х,0+^1(х,0, где ^1(х,0= =ук(х,0 для некоторого к, при этом индекс к «не фигурирует» в пакете Далее, пусть {х1, ..., xN.1} - набор узлов волны ук, для которых справедливо утверждение 1 по отношению к пакету Согласно индуктивному предположению (справедливость начала шага индукции согласно (10)) будет: Г1(^_1(хь-), ..., ^-1(х№1,0, dЪ1N-l(Xl,•)/dt, ..., d^N-1(xN-1,•)/dt)^0. Следовательно, имеет место: Г^^хл, •), ..., Ы%-1,-), ШЫтьОШ, ...,
Ш^(х№1,0Ш)^0, (12)
поскольку в противном случае найдутся вещественные числа г7,7=1,.,2(М-1), из которых не все равны нулю, такие, что (ниже 0 е С Ш(Т^)) выполнимо:
N-1(xN-1,•)/dt =
= Г 1^1 (х],•) -.... -.Г№1^1(х№1,0 - ^^(х^^Ш-....
-.r2(N-1)d^l(XN-1,•)/dt =
^-^1Ук(х1,•) -. ... -.г№1Ук(х№1,0 -.rNdуkСxl,•)/dt-...
-.r2(N-1)dуk(XN-1,•)/dt.
Данные равенства приводят к противоречию в силу уравнений (7) и положения, что индекс к не присутствует в (т. е. элементарные функции волны ук не входят в состав элементарных функций конструкции пакета ^л).
Далее, для стоячей волны через xN обозначим ее неузловую точку, не принадлежащую {х1, ..., xN-1}. Покажем, что пара (^№{хь ..., xN}) удовлетворяет утверждению 1. Будем рассуждать от противного. Пусть
Г^^ху), ..., Ых№0, d^N(Xl,•)/dt, ..., Ш^(х№0М)=0.
Тогда приходим к очевидной цепочке (ниже 0еС Ш(Т,Я)) отношений:
Г1^-1(хь-) +..
+ = Г1^1(х1,-) -...
Ук(х1,-)-.. .-^Ук^О-^+^ДХьОМ -...
где не все г7, 7=1,...ДК нулевые. Вариант г^0, г2М=0 исключаем согласно (12), случай г^0, г2М=0 противоречит соотношениям (7); вариант гм=0, г2^0 или г^0, г2^0 приводит к аналогичным рассуждениям. Наконец, вектор-функция СО1(^г(Х1,0, ..., Ы%,0, d^N(Хl,•)/d1, ..., d^N(xN,•)/d1),
1еТ
имеет дифференциальную реализацию (9) в силу (6), (7), (10).
З а м е ч а н и е 5: а) условие (11) не является характеристическим по отношению к утверждению 1 (поэтому его не использовали в доказательстве), но в Модуле-4 оно имеет важную роль при определении числа N
б) по существу свойство, сформулированное в утверждении 1, можно интерпретировать как структурный признак существования «квадратной» матрицы наблюдающего устройства у(1)=Сх(1), и такой, что при этом матрица С суть невырожденная (в конечном счете, в основе оттого положения лежит так называемое типическое свойство [19, с. 52].
4. Смежные вопросы
В заключение отметим, что естественным развитием программно-алгоритмической среды «МОДФОКБ» может служить задача имитационного моделирования форм-колебаний балки Тимошенко (Т1-балки; т. е. уточнение в прогибе балки неперпендикулярной деформации плоского сечения поперечного сдвига с учетом инерции его поворота [20]) с жестко закрепленным концом и присоединенным на свободном конце балки твердым телом. При этом в данном расширении программной среды «МОДФОКБ» функционально можно предусмотреть [21-24]:
1) задание физических параметров Т1-балки (модули упругости и сдвига, плотность материала балки, момент инерции и площадь сечения балки) и физических характеристик прикрепленного твердого тела (масса, момент инерции, статический момент инерции в торцевых координатах балки);
2) решение операторно-спектральной задачи (расчет собственных частот колебаний механической системы «Тьбалка + твердое тело»);
3) построение собственных изгибных форм (численное интегрирование дифференциальных уравнений (27) с начально-краевыми условиями (14)-(16));
4) адаптация параметров пунктов 1) - 3) к представлению в Модуль-1.
Заметим, что позиции 2) и 3) определяют решение задачи о собственных изгибных колебаниях Тьбалки с одним жестко закрепленным концом и с присоединенным на другом конце этой балки абсолютно твердым телом; для данной задачи существует [25] регулярный способ вычисления параметров собственных форм Тьбалки, основанный на редукции к задаче интегрирования обыкновенного дифференциального уравнения вида
- X2(1-X2r2q2)w(x)=0, (13)
при следующих граничных условиях: ^(0)=0, (q2diw(x)/dxъ+(1+qAX2)dw(x)/dx)x=0=0, (14) ((1 +q2LX2)d3w(x)/dx3+(( 1 +q4X2)L+ +(r2+q2))X2dw(x)/dx+(1-X2r2q2)X2mw(x))x=l=0, (15) (q2JX2dъw(x)/dxъ-(1-X2r2q2)d2w(x)/dx2+ +(1 +q4X2)JX2 dw(x)/dx+ +(1-X2r2q2)X2(L-q2)w(x))x=l=0; (16)
здесь L - статический момент инерции твердого тела относительно плоскости системы координат с началом в точке крепления тела к балке, J - момент инерции твердого тела, т - масса твердого тела, X - безразмерная собственная частота колебаний системы «Т1-балка + твердое тело», г2 -квадрат радиуса инерции балки (коэффициент, учитывающий инерцию поворота поперечного сечения балки), q2 - коэффициент, учитывающий влияние деформации сдвига.
З а м е ч а н и е 6. Если в (14)-(16) положить r2=q2=0, то, как частный случай, получим задачу о собственных колебаниях балки Эйлера - Бернули с присоединенным твердым телом, если пренебречь только членами, учитывающими деформации сдвига (т. е. q2=0), будем иметь балку Релея.
БИБЛИОГРАФИЧЕСКИЙ СПИСОК
1. Рубановский В. Н. Устойчивость установившихся движений сложных механических систем // Итоги науки и техники : общая механика. М. : ВИНИТИ, 1982. Т. 5. С. 62-134.
2. Гантмахер Ф. Р. Лекции по аналитической механике. М. : Физматгиз, 1960. 296 с.
3. Hablani H. B., Shrivastava S. K. Nontrivial Equilibrium and Lyapunov Stability of Flexible Damped Gravity Satellite as a Hybrid Dynamical System // AIAA Pap. 1978. № 53. P. 1-10.
4. Сергиенко И. В., Дейнека В. С. Идентификация параметров динамической задачи теории упругости тела с включением // Кибернетика и системный анализ. 2009. № 3. С. 75-97.
5. Сомов Е. И., Бутырин С. А. Полетная геометрическая идентификация и калибровка космического телескопа и системы звездных датчиков // Идентификация систем и задачи управления тр. VIII Междунар. конф. / ИПУ РАН. М. 2009. С. 1189-1201. (Труды SICPRO'09).
6. Лурье А. И. Аналитическая механика. М. : Физматгиз, 1961. 824 с.
7. Румянцев В. В. Некоторые задачи динамики сложных систем // Проблемы прикладной математики и механики. М., 1971. С. 179-188.
8. Владимиров В. С. Уравнения математической физики. М. : Наука, 1981. 512 с.
9. Русанов В. А., Шарпинский Д. Ю. К теории структурной идентификации нелинейных многомерных систем // Прикладная математика и механика. 2010. Т. 74, вып. 1. С. 119-132.
10. Калман Р., Фалб П., Арбиб М. Очерки по математической теории систем. М. : Мир, 1971. 400 с.
11. Механика больших космических конструкций / Ба-ничук Н. В., Карпов И. И., Климов Д. М. [и др.]. М. : Факториал, 1997. 302 с.
12. Кузьмин П. А. Малые колебания и устойчивость движения. М. : Наука, 1973. 208 с.
13. Годунов С. К. Уравнения математической физики. М. : Наука, 1971. 416 с.
14. Гантмахер Ф. Р. Теория матриц. М. : Наука, 1988. 550 с.
15. Русанов В. А., Шарпинский Д. Ю., Козырев В. А., Удилов Т. В. Реализация динамической модели «РЕДИМ» // Свидетельство Федеральной службы по интеллектуальной собственности и патентам о гос. регистрации программы для ЭВМ № 2008612596. 2008.
16. Полетыкин А. Г., Байбулатов А. А. Обзор зарубежного опыта разработки систем человеко-машинного интерфейса АСУ ТП // Идентификация систем и задачи управления : тр. III Междунар. конф. / ИПУ РАН. 2004. С. 719-734. - (Труды SICPRO'04).
17. Колмогоров А. Н., Фомин С. В. Элементы теории функций и функционального анализа. М. : Наука, 1976. 544 с.
18. Данеев А. В., Русанов В. А. Геометрические характеристики свойств существования конечномерных (А,В)-моделей в задачах структурно -параметрической идентификации // Автоматика и телемеханика. 1999. № 1. С. 3-8.
19. Уонэм М. Линейные многомерные системы управления : Геометрический подход. М. : Наука, 1980. 376 с.
20. Тимошенко С. П., Гудъер Дж. Теория упругости. М. : Наука, 1979. 560 с.
21. Парафесь С. Г. Метод идентификации тонкостенных конструкций // Идентификация систем и задачи управления : тр. III Междунар. конф. / ИПУ РАН. 2004. С. 326-335. (Труды SICPRO'04).
22. Калинин Е. Н. Идентификация модели упругой роторной системы // Идентификация систем и задачи управления : тр. III Междунар. конф. / ИПУ РАН. 2004. С. 653-659. (Труды SICPRO'04).
23. Моделирование распределённых систем реального времени / Галатенко В. А., Костюхин К. А., Шмырев Н. В., Малиновский А. С. // Параллельные вычисления и задачи управления : тр. IV Междунар. конф. / ИПУ РАН. 2008. С. 1349-1355. (Труды РАС0'08).
24. Копысов О. Ю., Кулагин В. П. К проблеме идентификации параметров динамики вращения твердого тела // Идентификация систем и задачи управления : тр. VIII Междунар. конф. / ИПУ РАН. 2009. С. 492500. (Труды SICPRO'09).
25. Ахтямов А. М., Урманчеев С. Ф. Определение параметров твердого тела, прикрепленного к одному из концов балки, по собственным частотам колебаний // Сибирский журн. индустриал. математики. 2008. Т. XI, № 4. С. 19-24.
