Прикладная математика и механика, 2023, T. 87, № 5, стр. 801-819

1. Введение. Среди моделей, предложенных для приближенного решения начально-краевых задач механики упругих протяженных тел (стержней, валов и т.п.), особое место занимает теория балок, основанная на выдвинутых Я. Бернулли гипотезах [1] о том, что во время движения поперечное сечение не деформируется и остается ортогональным центральной линии тела. Хотя эта теория применима для широкого класса задач, она не учитывает влияние на состояние деформируемого тела сдвиговых перемещений, депланации и деформации поперечного сечения, связи продольных и поперечных движений, обусловленных эффектом Пуассона, и т.п. Для компенсации этих дефектов модели применяются уточняющие формулы, позволяющие учесть сжатие или растяжение поперечных сечений при продольных перемещениях (поправка Релея [2]), а также сдвиги и инерционность поворотов сечения при упругом изгибе (модель балки Тимошенко [3]). В модели Райсснера [4] применен вариационный подход для вывода уравнений, описывающих упругий изгиб тонкой пластины (балки), для определенного заранее распределения полей перемещений в поперечном направлении. Вариационные формулировки используются для вывода совместных балочных уравнений и более высокого порядка, в которых различным образом можно учесть пространственное распределение как перемещений, так и напряжений в упругом теле [5].

Цель этой работы – развитие метода интегродифференциальных соотношений (МИДС) [6], который основан на интегральном представлении законов состояния механической системы, в применении к задачам управления упругими слоистыми системами. Основным преимуществом МИДС является то, что поля перемещений и напряжений, которые являются независимыми переменными, строго удовлетворяют уравнению динамического равновесия и граничным условиям. В явном виде могут быть получены двусторонние энергетические оценки качества любого допустимого движения. Ранее на основе МИДС разработаны регулярные проекционные и вариационные подходы к построению математических моделей, описывающих движения упругих стержней различного сечения [7, 8]. Для расчета спектральных характеристик неоднородных валов на основе МИДС [9] применялся метод ускоренной сходимости [10].

Управление колебаниями протяженных упругих (стержневых) систем можно осуществлять как посредством приложения сил или моментов на концевых торцах, так и путем нагрузки по длине. Последняя может создаваться, например, пьезоэлектрическими актюаторами (элементами). Зачастую предполагается для простоты, что актюаторы размещены по длине стержня периодически, хотя возможно рассмотреть и более общие геометрии [11], а также задачи оптимизации расположения актюаторов [12, 13]. В работах [14, 15] затронут вопрос управляемости отдельных мод продольных колебаний упругих стержней и разбиения мод на группы, каждая из которых управляется определенной комбинацией входов. Распространенный подход при моделировании стержневых систем, управляемых пьезоэлементами, состоит в добавлении к уравнениям изгиба правых частей, содержащих дельта-функции [16], отражающие наличие скачков напряжений в стержне. Однако такой подход может приводить к трудностям в моделировании вследствие необходимости учитывать эти обобщенные функции в расчетах, а также может понизить класс решения.

В представленной работе применяется новый подход, состоящий в постановке вариационной задачи для описания напряженно-деформированного состояния как в балке, так и в актюаторах. Поскольку сами пьезоэлементы представляют собой упругие тела, представляется оправданным моделирование системы “балка-пьезоэлемент” как слоистое тело [17, 18], в котором упругие напряжения в пьзоэлектрических слоях заданы как функции времени (например, вследствие некоторого подаваемого электрического напряжения), а поля напряжений и перемещений внутри балки генерируются вследствие подходящих краевых и межэлементных условий. На основе МИДС для расчета трехмерной динамики системы “балка-пьезоэлементы”, который состоит из центрального упругого слоя и двух граничных пьезоэлектрических слоев, разбитых на последовательность пар работающих антисимметрично актюаторов, выбрана полиномиальная по поперечным координатам аппроксимация перемещений и напряжений.

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

2. Описание динамической системы. Рассматриваются изгибные движения тонкого прямолинейного упругого стержня с прямоугольным сечением. Декартова система координат $O{{x}_{1}}{{x}_{2}}{{x}_{3}}$ выбирается так, что начало координат находится в центре одного из торцевых сечений, ось ${{x}_{1}}$ направлена в центр противолежащего торца, оси ${{x}_{2}}$ и ${{x}_{3}}$ параллельны сторонам поперечного сечения. Обозначим пространственную область, занимаемую стержнем в недеформированном состоянии, как

где $a$, $2b \ll a$, $2c \ll a$ – линейные размеры тела.

К боковым поверхностям ${{x}_{2}} = \pm b$ попарно без зазоров присоединены пьезоэлектрические актюаторы (элементы). Схематически разрез стержня с прикрепленными элементами в плоскости $O{{x}_{1}}{{x}_{2}}$ изображен на рис. 1. В недеформированном состоянии $2N$ элементов занимают соответствующие пространственные области:

где ${{\chi }_{n}}$ – продольные координаты границ пьезоэлементов. В каждой паре пьезоэлектрических элементов с помощью антисимметричного электрического потенциала создаются однородные механические напряжения $ \pm {{s}_{n}}(t)$, управляющие на временном промежутке $t \in \mathcal{T} = (0,T)$ изгибом стержня (см. рис. 2). Введем соответствующее тензорное поле пьезоэлектрических напряжений

Здесь и далее ${{{\mathbf{e}}}_{j}}$ – орт вдоль оси ${{x}_{j}}$, $j = 1,2,3$.

Обозначим полный объем, занимаемый этой электромеханической системой в пространстве, и соответствующую пространственно-временную область как

Положим для простоты, что в отсутствие электрического потенциала как материал стержня, так и материал пьезоактюаторов однороден и изотропен. Используются следующие сравнимые для обоих материалов константы: объемная плотность ${{\rho }_{j}}$, модуль Юнга ${{E}_{j}}$, коэффициент Пуассона ${{\nu }_{j}}$, $j = 0,1$. Материальные константы стержня и управляющих элементов помечены соответственно нижними индексами 0 и 1. Симметрический тензор модулей податливости четвертого ранга ${\mathbf{S}}({\mathbf{x}})$ и функция плотности $\rho ({\mathbf{x}})$ полностью задаются на $\mathcal{X}$ этими шестью параметрами.

В предлагаемой слабой формулировке начально-краевой задачи линейной теории упругости искомыми функциями являются вектор перемещений ${\mathbf{v}}:\mathcal{D} \to {{\mathbb{C}}^{3}}$ и симметрический тензор второго ранга $\eta :\mathcal{D} \to {{\mathbb{C}}^{{3 \times 3}}}$, $\eta = {{\eta }^{{\text{T}}}}$, который представляет интеграл по $t$ от тензора механических напряжений. В выбранной постановке удобно использовать комплексные переменные для нахождения собственных значений и векторов задачи, но физический смысл имеют лишь действительные части вектора ${\mathbf{v}}({\mathbf{x}},t)$ и тензора $\eta ({\mathbf{x}},t)$. Для сокращения записи неизвестные вектор и тензор объединяются в одну функцию ${\mathbf{V}} = ({\mathbf{v}},\eta )$, которая предполагается квадратично интегрируемой со всеми своими первыми производными: ${\mathbf{V}} \in {{H}^{1}}(\mathcal{D};{{\mathbb{C}}^{9}})$.

Обобщенная частная производная по времени тензора $\eta ({\mathbf{x}},t)$ определяет напряжения $\sigma ({\mathbf{x}},t)$, возникающие в материале, а дивергенция – объемную плотность импульса ${\mathbf{p}}({\mathbf{x}},t)$ в форме

Для таким образом введенных переменных как перемещения $v$, так и напряжения $\sigma $ – соболевские функции по пространственным переменным.

С учетом введенных далее граничных условий после подстановки (2.5) автоматически (в обобщенной форме) выполняется второй закон Ньютона:

В уравнении динамики (2.6) введена тестовая вектор-функция ${\mathbf{w}}:\mathcal{D} \to {{\mathbb{C}}^{3}}$, которая вместе с тензорным полем $\zeta :\mathcal{D} \to {{\mathbb{C}}^{{3 \times 3}}}$, $\zeta = {{\zeta }^{{\text{T}}}}$, составляют используемую в дальнейшем тестовую функцию ${\mathbf{W}} = ({\mathbf{w}},\zeta )$.

Используя переменную ${\mathbf{V}}$, классические уравнения состояния линейной упругости

можно переписать в интегральной форме: где введен линейный оператор и билинейная форма:

В качестве граничных условий на поля ${\mathbf{V}}$ и ${\mathbf{W}}$ выбираются нулевые напряжения на боковых сторонах стержня и шарнирное закрепление торцов:

Здесь интегралы по поперечному сечению $\mathcal{S}$ означают момент сил относительно оси ${{x}_{3}}$ в сечении с координатой ${{x}_{1}}$, а ${{\alpha }_{1}}$, ${{\beta }_{1}}$ – множители Лагранжа, являющиеся произвольными комплексными функциями времени, которые означают углы поворота торцевых сечений при ${{x}_{1}} = 0$ и ${{x}_{1}} = a$. Пространство $H$ в (2.8) определяется как пространство соболевских функций, чьи следы удовлетворяют (2.10).

Известные начальные распределения перемещений ${\mathbf{v}}$ и плотности импульса ${\mathbf{p}}$ задаются соотношениями:

3. Задача на собственные значения. Начально-краевой задаче (2.8)–(2.11) соответствует краевая задача на собственные частоты $\omega \in \mathbb{R}$:

с граничными условиями где ${{\tilde {\alpha }}_{1}}(0)$, ${{\tilde {\alpha }}_{1}}(a)$, ${{\tilde {\beta }}_{1}}(0)$, ${{\tilde {\beta }}_{1}}(a)$ – произвольные комплексные константы, а $\tilde {H}$ – подпространство соболевского пространства с ограничениями (3.2).

4. Аппроксимация полей перемещений и напряжений. Для нахождения собственных частот и форм колебаний используем полиномиальные аппроксимации векторных и тензорных полей ${\mathbf{\tilde {v}}}$, ${\mathbf{\tilde {w}}}$, $\tilde {\eta }$, $\tilde {\zeta }$ относительно поперечных координат ${{x}_{2}}$, ${{x}_{3}}$:

Здесь ${{{\mathbf{v}}}^{{(jk)}}}({{x}_{1}})$, ${{\eta }^{{(jk)}}}({{x}_{1}})$ – неизвестные векторно- и тензорнозначные пробные функции соответствующей размерности, ${{{\mathbf{w}}}^{{(jk)}}}({{x}_{1}})$, ${{\zeta }^{{(jk)}}}({{x}_{1}})$ – произвольные тестовые функции. Заметим, что для более точного учета слоистой структуры рассматриваемой системы можно было бы использовать отдельные разложения в каждом из трех слоев (${{x}_{2}} \in ( - h, - b)$, ${{x}_{2}} \in ( - b,b)$, ${{x}_{2}} \in (b,h)$), учитывая соответствующие условия склейки перемещений и напряжений при ${{x}_{2}} = \pm b$, но для построения простейшей модели изгиба ограничимся общей аппроксимацией по всему поперечному сечению $\mathcal{S}$.

В [7] показано, что для стержней имеющих две плоскости геометрической и механической симметрии существует по крайней мере четыре независимых типа движения: кручение, продольные колебания, а также два изгиба относительно осей ${{x}_{2}}$ и ${{x}_{3}}$. Как указано в табл. 1, они отличаются определенной четностью или нечетностью полиномов по ${{x}_{2}}$ и ${{x}_{3}}$ в (4.1). Кручение выделяется тем, что точки центральной линии стержня (${{x}_{2}} = {{x}_{3}} = 0$) всегда остаются неподвижными, а сечение поворачивается вокруг оси ${{x}_{1}}$. Продольные движения оставляют точки центральной линии на оси $O{{x}_{1}}$. При изгибе вокруг оси ${{x}_{2}}$ центральная линия смещается вдоль оси ${{x}_{3}}$. И, наоборот, при изгибе вокруг ${{x}_{3}}$ центральная линия перемещается только вдоль ${{x}_{2}}$.

Так как тензор пьезоэлектрических напряжений $\xi $ антисимметричен относительно плоскости $O{{x}_{1}}{{x}_{3}}$ и симметричен относительно $O{{x}_{1}}{{x}_{2}}$, он возбуждает только изгибные деформации вокруг оси ${{x}_{3}}$. Ограничимся рассмотрением лишь этого типа колебаний и построим приближенные поля перемещений и напряжений ограничив максимальные степени разложения по четным или нечетным степеням ${{x}_{2}}$ и ${{x}_{3}}$. Для получения совместной краевой задачи на собственные значения относительно ${{x}_{1}}$ необходимо правильно подобрать как подходящие степени, так и максимальный порядок разложения каждой из компонент ${\mathbf{\tilde {v}}}$ и $\tilde {\eta }$. В результате получаем допустимую аппроксимацию:

где ${{N}_{p}}$ – порядок приближения. Эти аппроксимации автоматически выполняют граничные условия (3.2) при ${{x}_{2}} = \pm h$ и ${{x}_{3}} = \pm c$. Приближения для ${\mathbf{\tilde {w}}}$ и $\tilde {\zeta }$ аналогичны (4.2) и для сокращения объема не приводятся.

Далее ограничимся полиномиальными разложениями компонент полей ${\mathbf{\tilde {v}}}$, ${\mathbf{\tilde {w}}}$, $\tilde {\eta }$, $\tilde {\zeta }$ наименьшей степени: ${{N}_{p}} = 0$. Учитывая, что индексы $j$, $k$ в (4.2) принимают только неотрицательные значения, и вводя новые обозначения для неизвестных функций от ${{x}_{1}}$, получаем

Граничные условия на торцах ${{x}_{1}} = 0$ и ${{x}_{1}} = a$ переходят в смешанные краевые ограничения

Подставляя эти вектор-функции в интегральное соотношение (3.1) и учитывая произвольность выбора ${\mathbf{r}}$, относительно ${{x}_{1}}$ получим систему линейных дифференциально-алгебраических уравнений четвертого порядка с действительными коэффициентами:

Искомые пробные функции ${{q}_{j}}({{x}_{1}}) \in \mathbb{R}$ и произвольные тестовые функции ${{r}_{j}}({{x}_{1}}) \in \mathbb{R}$ собираются соответственно в векторы ${\mathbf{q}}({{x}_{1}}) \in {{\mathbb{R}}^{6}}$ и ${\mathbf{r}}({{x}_{1}}) \in {{\mathbb{R}}^{6}}$. В векторном виде это уравнение можно представить как

Компоненты постоянных действительных матриц ${\mathbf{B}} \in {{\mathbb{R}}^{{6 \times 6}}}$ и ${\mathbf{C}} \in {{\mathbb{R}}^{{6 \times 6}}}$ определяются аналитически и зависят от геометрических и механических параметров упругой системы $b$, $c$, $h$, ${{\rho }_{j}}$, ${{E}_{j}}$, ${{\nu }_{j}}$, $j = 0,1$.

Последние два алгебраических уравнения в (4.5) всегда можно разрешить относительно ${{q}_{5}}$ и ${{q}_{6}}$. После их исключения, характеристическое уравнение первых четырех дифференциальных уравнений в (4.5) уже только относительно кинематических переменных ${{q}_{1}}$, ${{q}_{2}}$ и динамических переменных ${{q}_{3}}$, ${{q}_{4}}$ имеет вид

где ${{d}_{j}} > 0$, $j = 1,2,3$, постоянные коэффициенты, вычисляемые алгебраически через компоненты матриц ${\mathbf{B}}$ и ${\mathbf{C}}$.

Учитывая, что частота $\omega $ принимает только действительные значения, можно, разрешив биквадратное уравнение (4.7), выразить в явном виде волновое число $\kappa $ через квадрат частоты:

Анализ показывает, что вне зависимости от параметров системы выполняется строгое неравенство $d_{1}^{2} > {{d}_{3}}$ и квадрат волнового числа ${{\kappa }^{2}}({{\omega }^{2}})$ в (4.8) принимает только действительные значения. Один из корней всегда отрицателен, второй тоже больше нуля, если ${{\omega }^{2}} < \omega _{*}^{2} = {{d}_{2}}{\text{/}}{{d}_{3}}$. Для высоких частот (${{\omega }^{2}} > \omega _{*}^{2}$) оба значения отрицательны (${{\kappa }^{2}} < 0$), а, следовательно, и все четыре корня характеристического уравнения (4.7) чисто мнимые.

Будем искать решение краевой задачи (4.4)–(4.5) в виде

Можно непосредственно убедиться, что компоненты собственного вектора ${{{\mathbf{q}}}_{m}}({{x}_{1}})$ удовлетворяют краевым условиям (4.4). Все параметры $q_{m}^{j}$, ($j = 1, \ldots ,6$) этого вектора с точностью до множителя явно находятся как алгебраические функции от собственной частоты ${{\omega }_{m}}$ после подстановки ${{{\mathbf{q}}}_{m}}$ и ${{\omega }_{m}}$ в (4.5). Волновое число ${{\kappa }_{m}} > 0$ и частота ${{\omega }_{m}} \in \mathbb{R}$ удовлетворяют характеристическое уравнение (4.7) после подстановки в него чисто мнимого числа $\kappa = \pm i{{\kappa }_{m}}$ и частоты $\omega = {{\omega }_{m}}$.

Так как $d_{1}^{2} > {{d}_{3}}$, для каждого натурального числа $m$ существуют два типа собственных движений упругой системы, которые определяются выбором плюса (верхняя ветвь колебаний) или минуса (нижняя ветвь) в выражении для $\omega _{m}^{2}$ из (4.9). Положим в дальнейшем, что ${{\omega }_{{ \pm m}}} > 0$ положительные значения частот для верхней и нижней ветви соответственно. Как собственное значение подходят также и отрицательные частоты $ - {{\omega }_{{ \pm m}}}$. В результате одно натуральное число $m$ порождает четыре собственных вектора

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

Для выбора неопределенного параметра собственного вектора ${{{\mathbf{q}}}_{{ \pm m}}}$ естественно зафиксировать его энергетическую норму. Согласно (4.3), функция ${\mathbf{q}}_{{ \pm m}}^{ + }({{x}_{1}})$ определяет в пространственной области $\mathcal{X}$ векторное кинематическое поле ${{{\mathbf{v}}}_{{ \pm m}}}$ и тензорное динамическое поле ${{\eta }_{{ \pm m}}}$:

Так как представленные поля могут принимать комплексные значения, для получения действительного решения начально-краевой задачи (2.8)–(2.11) нам потребуется как собственный вектор ${{{\mathbf{V}}}_{{ \pm m}}}({\mathbf{x}})$, соответствующий положительной частоте ${{\omega }_{{ \pm m}}} > 0$ и вектору ${\mathbf{q}}_{{ \pm m}}^{ + }$, так и сопряженный ему вектор ${{\overline {\mathbf{V}} }_{{ \pm m}}}({\mathbf{x}})$, соответствующий отрицательной частоте $ - {{\omega }_{{ \pm m}}}$ и вектору ${\mathbf{q}}_{{ \pm m}}^{ - }$.

Нормировочное соотношение для полной механической энергии моды с индексом $ \pm m$ имеет вид

Отметим, что кинетическая и потенциальная энергия собственных колебаний равны. Знак сопряжения над ${{\eta }_{{ \pm m}}}$ в выражении для плотности потенциальной энергии выбран из-за того, что собственный тензор ${{\eta }_{{ \pm m}}}$, согласно (4.3), принимает мнимые значения. Важно подчеркнуть также, что различные собственные векторы ортогональны относительно скалярного произведения, порождающего энергетическую норму:

5. Пример расчета собственных частот и форм колебаний. В качестве примера рассмотрим колебания упругой системы со следующими безразмерными геометрическими и механическими параметрами

Подставим эти параметры вместе с соответствующими им аппроксимациями (4.3) в интегральное уравнение (3.1). После интегрирования по координатам ${{x}_{2}}$ и ${{x}_{3}}$ находятся компоненты матриц ${\mathbf{B}}$ и ${\mathbf{C}}$, а из них – коэффициенты характеристического уравнения (4.7).

Построим амплитудно-частотные характеристики системы с выбранными параметрами (5.1). На рис. 3 изображен график зависимости абсолютных значений корней ${\text{|}}\kappa {\text{|}}$ уравнения (4.7) от частоты $\omega > 0$. Четыре корня можно представить либо как $\kappa \in \{ i{{\kappa }^{0}}, - i{{\kappa }^{0}},{{\kappa }^{1}}, - {{\kappa }^{1}}\} $ при $\omega < {{\omega }_{*}}$, либо как $\kappa \in \{ i{{\kappa }^{0}}, - i{{\kappa }^{0}},i{{\kappa }^{1}}, - i{{\kappa }^{1}}\} $ при $\omega > {{\omega }_{*}}$. На графике критическая частота ${{\omega }_{*}}$ помечена вертикальной пунктирной линией и в нашем примере равна:

Действительные корни характеристического уравнения $ \pm {{\kappa }^{1}}(\omega )$ в докритической зоне $\omega \in (0,{{\omega }_{*}})$ отображаются на рис. 3 штриховой кривой и для шарнирного закрепления концов не реализуются в общем решении (4.9). На этом интервале функция ${{\kappa }^{1}}(\omega )$ достигает максимума ${{\kappa }_{1}} = \kappa {\kern 1pt} * \approx 0.778$ при $\omega {\kern 1pt} * \approx 0.467$. В отличие от модели балки Бернулли, в рамках которой абсолютные значения всех корней совпадают, в предложенной модели всегда ${{\kappa }^{1}}(\omega ) < {{\kappa }^{0}}(\omega )$ для допустимого интервала $(0,{{\omega }_{*}})$.

Для сравнения штрихпунктирной линией на рисунке построен график квадратичной зависимости от частоты волнового числа ${{\kappa }^{c}}(\omega ) = \kappa _{0}^{c}{{\omega }^{2}}$, $\kappa _{0}^{c} \approx 0.939$, классической модели балки c соответствующими осредненными по сечению параметрами. При $\omega \to 0$ три функции ${{\kappa }^{b}}$, ${{\kappa }^{0}}$, ${{\kappa }^{1}}$ (штрихпунктирная, сплошная и штриховая линии вблизи начала координат на рис. 3) асимптотически склеиваются. При выбранном отношении высоты стержня $2h$ к его длине $a$ ($1:10$), графики заметно расходятся уже при достаточно низких частотах ($\omega \sim 0.3$).

Если $\omega = {{\omega }_{*}}$, то ${{\kappa }^{1}} = 0$, а в закритической области $\omega \in ({{\omega }_{*}},\infty )$ появляются два мнимых корня, абсолютное значение которых ${{\kappa }^{1}}(\omega )$ показано в нижней правой части графика сплошной линией. Это строго монотонно возрастающая функция, которая всегда ниже абсолютных значений двух других сопряженных корней (${{\kappa }^{1}}(\omega ) < {{\kappa }^{0}}(\omega )$). Ветвь ${{\kappa }^{0}}(\omega )$ также строго монотонна и изображена сплошной линией в верхней левой части рис. 3.

Амплитудно-частотные характеристики, описанные выше, верны для любых консервативных типов закрепления стержня на торцах. При выбранной шарнирной заделке обоих концов собственные частоты и векторы находятся в виде (4.9). Для построения таких форм колебаний используются только мнимые волновые числа. Выбирая полученные частоты ${{\omega }_{m}}$ и волновые числа ${{\kappa }_{m}}$, отобразим их на рис. 3 как точки плоскости $(\omega ,\kappa )$. Кружками обозначена нижняя ветвь колебаний $({{\omega }_{{ - m}}},{{\kappa }_{m}})$, а ромбами – верхняя $({{\omega }_{{ + m}}},{{\kappa }_{m}})$. Как видно на графике эти точки строго ложатся на мнимые ветви функций ${{\kappa }^{0}}(\omega )$, $\omega > 0$, и ${{\kappa }^{1}}(\omega )$, $\omega > {{\omega }_{*}}$.

Для более точного частотного анализа можно обратиться к табл. 2, в которую помещены значения первых десяти частот $\omega _{m}^{c}$, $m = 1, \ldots ,10$, изгиба шарнирно закрепленной балки Бернулли с параметрами соответствующими рассматриваемой в этом примере упругой системы (вторая строка). Третья строка табл. 2 соответствует нижней ветви собственных частот ${{\omega }_{{ - m}}}$, которые получены в рамках предложенной модели. В последней строке указаны десять собственных частот ${{\omega }_{{ + m}}}$ верхней ветви.

Небольшое расхождение $\omega _{1}^{c}$ и ${{\omega }_{{ - 1}}}$ (порядка 3%) видно даже для первой моды колебаний. С возрастанием номера $m$ это рассогласование быстро нарастает достигая почти 70% для $m = 10$. При этом классическая частота $\omega _{{10}}^{c}$ превышает критическую ${{\omega }_{*}}$, выше которой появляются собственные частоты верхней ветви ${{\omega }_{{ + m}}}$. Их величины на два порядка выше основной частоты ${{\omega }_{{ - 1}}}$, но в относительном выражении растут гораздо медленнее с номером $m$ чем ${{\omega }_{{ - m}}}$.

Рассмотрим соответствующие двум ветвям $({{\omega }_{{ \pm m}}},{{\kappa }_{m}})$ формы собственных колебаний. На рис. 4 показаны углы упругих поворотов в заделке ${{x}_{1}} = 0$ для различных $m$:

где ${{{v}}_{{j, \pm m}}}$, $j = 1,2,3$, – компоненты вектора собственных перемещений ${{{\mathbf{v}}}_{{ \pm m}}}$. Значение ${{\alpha }_{1}}$ – отмеряемый от недеформированного положения угол поворота поперечного сечения, ${{\alpha }_{2}}$ – угол наклона деформированной центральной линии стержня (см. рис. 2). Оба поворота осуществляются вокруг оси ${{x}_{3}}$. На рис. 4 значения ${{\alpha }_{1}}( - m)$ для нижней ветви колебаний отмечены закрашенными ромбами, ${{\alpha }_{1}}( + m)$ для верхней – незакрашенными. Углы ${{\alpha }_{2}}( - m)$ для нижней ветви показаны закрашенными кружками, ${{\alpha }_{2}}( + m)$ – незакрашенными. Не должны смущать большие значения этих параметров, которые соответствуют собственным колебаниям с безразмерной энергией, равной единице. Это связано с тем что модуль Юнга материала стержня равен единице, а у пьезоматериала этот модуль еще меньше. Для анализа собственных форм важны только относительные значения этих углов.

Как видно из графика, для нижней ветви колебаний, как угол поворота сечения ${{\alpha }_{1}}$, так и угол наклона центральной линии ${{\alpha }_{2}}$ положительны и монотонно уменьшаются с увеличением номера моды $m$. Разница между этими углами равна углу сдвига ${{\gamma }_{{ - m}}} = {{\alpha }_{2}}( - m) - {{\alpha }_{1}}( - m)$ в центре торцевого сечения; он увеличивается для больших $m$, что указывает на усиление роли перерезывающих сил ${{\sigma }_{{12, - m}}}$ в этом типе колебаний. Для верхней ветви ($ + m$) важно отметить, что ${{\alpha }_{1}}$ и ${{\alpha }_{2}}$ разных знаков и почти линейно растут. При малых номерах $m$ значения ${{\alpha }_{2}}( + m)$ близки к нулю по сравнению с ${{\alpha }_{2}}( - m)$, а сдвиговые деформации ${{\gamma }_{{ - m}}}$ – существенны, что выявляет сдвиговый характер этих колебаний в отличие от изгибного типа для нижней ветви.

Двумя другими параметрами собственных движений, определяемыми из дифференциальных соотношений, являются амплитуды нормальных и сдвиговых напряжений ${{\sigma }_{{11}}}$ и ${{\sigma }_{{12}}}$. На рис. 5 показаны максимальные значения нормальных и сдвиговых напряжения для различных $m$:

При этом напряжения ${{\sigma }_{1}}$ помечены ромбами, ${{\sigma }_{2}}$ – кружками; нижняя ветвь обозначена закрашенными символами, верхняя – незакрашенными. Для нижней ветви с ростом номера $m$ напряжение ${{\sigma }_{1}}$ уменьшается, ${{\sigma }_{2}}$ увеличивается, при этом ${{\sigma }_{1}} > {{\sigma }_{2}}$. Для низших частот этой ветви амплитуда нормальных напряжений в сечении ${{\sigma }_{1}}$ много больше сдвиговых ${{\sigma }_{2}}$, что указывает на преимущественно изгибный характер колебаний. Для верхней же ветви ${{\sigma }_{1}}( + m)$ – монотонно возрастающая функция, а ${{\sigma }_{2}}( + m)$ –убывающая. При малых $m$ сдвиговые напряжения много больше по амплитуде чем нормальные, это означает существенно сдвиговый характер собственных напряжений.

6. Декомпозиция начально-краевой задачи. Согласно методу Фурье представим пробные и тестовые переменные ${\mathbf{V}}$, ${\mathbf{W}}$ начально-краевой задачи (2.8)–(2.11) как разложение по собственным векторам задачи (4.4), (4.6) на аппроксимациях (4.3):

Здесь ${{y}_{j}}$, ${{z}_{j}}$ – неизвестные функции времени, $\delta {{y}_{j}}$, $\delta {{z}_{j}}$ – их вариации.

Подставим выражения для ${\mathbf{V}}$, ${\mathbf{W}}$ из (6.1) в (2.7). Учитывая условия (4.12), (4.13) для ортонормированного базиса ${{\{ {{{\mathbf{V}}}_{j}},{{\overline {\mathbf{V}} }_{j}}\} }_{{j \in \mathbb{Z}\backslash \{ 0\} }}}$, проинтегрируем выражения под пространственно-временными интегралами в (2.9) только по объему $\mathcal{X}$. Помня о произвольности вариаций $\delta {{y}_{j}}$, $\delta {{z}_{j}}$, получаем счетную систему обыкновенных уравнений первого порядка:

где индексу $j = - m$, $m \in \mathbb{N}$, соответствует собственная частота нижней ветви колебаний ${{\omega }_{{ - m}}} > 0$, а $j = + m$ – частота верхней ветви ${{\omega }_{{ + m}}} > {{\omega }_{{ - m}}}$.

Однородные граничные условия (2.10) для приближения ${\mathbf{V}}({\mathbf{x}},t)$, выбранного в виде (6.1), выполняются точно, так как собственные векторы ${{{\mathbf{V}}}_{j}}({\mathbf{x}})$ не нарушают ограничения краевой задачи (3.2). Чтобы получить начальные условия для неизвестных ${{y}_{j}}$, ${{z}_{j}}$ следует спроецировать заданные в момент $t = 0$ распределения ${{{\mathbf{v}}}^{0}}({\mathbf{x}})$ и ${{{\mathbf{p}}}^{0}}({\mathbf{x}})$ из (2.11) соответственно на собственные функции ${{{\mathbf{v}}}_{j}}({\mathbf{x}})$ и $\nabla \cdot {{\eta }_{j}}({\mathbf{x}})$. В результате получим систему начальных условий:

Из структуры задачи Коши (6.2), (6.3) сразу следует, что ${{z}_{j}}(t) = {{\bar {y}}_{j}}(t)$ и приближенное решение в виде (6.1) принимает только действительное значение: $\Im {\mathbf{V}}({\mathbf{x}},t) = 0$.

Правые части уравнений (6.2) имеют важную особенность для выбранного (шарнирного) способа закрепления торцов стержня. Напомним, что тензор пьезоэлектрических напряжений $\xi ({\mathbf{v}},t)$ имеет в заданной системе координат только одну ненулевую компоненту ${{\xi }_{{11}}}$, которая однородно распределена и независимо изменяется как функция времени ${{\xi }_{{11}}}({\mathbf{x}},t) = \pm {{s}_{n}}(t)$ на каждом пьезоэлементе $\mathcal{X}_{n}^{ \pm }$ и равна нулю на ${{\mathcal{X}}_{0}}$. Следовательно, правая часть, например, первого уравнения в (6.2) для индекса $j = \pm m$ вычисляется как

где координаты границ пьезоэлементов ${{\chi }_{n}}$ введены в (2.2), множитель ${{g}_{{ \pm m}}} \in \mathbb{R}$, общий для всех $n$, получается после интегрирования диагональных компонент ${{\eta }_{{kk, \pm m}}}({\mathbf{x}})$, $k = 1,2,3$, собственного тензора ${{\eta }_{{ \pm m}}}$ по части поперечного сечения

Из-за антисимметрии этих компонент напряжений возникает двойка перед знаком интеграла в первой строке (6.4) вместо дополнительного интегрирования по области $\mathcal{X}_{n}^{ - }$.

Покажем, что правая часть уравнений (6.2) – это линейная комбинация $N$ независимых функций времени с коэффициентами

которые принимают значения из конечного набора. Пусть номер одной моды равен ${{m}_{1}} = k$, $k = 0, \ldots ,N$, а другой ${{m}_{2}} = 2Nj \pm k$, $j \in \mathbb{N}$. Тогда для каждого номера пьезоэлемента $n$ выражение в квадратных скобках в (6.5) не меняется при переходе от ${{m}_{1}}$ к ${{m}_{2}}$, то есть правые части уравнений (6.2) с индексами $ \pm {{m}_{1}}$ и $ \pm {{m}_{2}}$ линейно зависимы. Важно отметить, что нулевая мода для шарнирного закрепления – это тривиальное решение, но соответствующая линейная комбинация тоже тривиальна и для каждой моды с номером $2Nj$ правая часть равна нулю и эти моды не управляемы. Итак, колебательные моды системы разбиваются на $N + 1$ классов эквивалентности по признаку линейной зависимости правых частей.

Введем новые функции управления

После такой замены задача Коши (6.2), (6.3) для определения неизвестной функции приобретает вид

7. Пример помодового управления движением стержня. Можно предложить различные стратегии управления системой (6.7), учитывая ее декомпозицию на группы мод. Исключая неуправляемую нулевую группу, каждая из остальных $N$ групп представляют собой счетное множество линейных осцилляторов с одним скалярным управлением ${{u}_{k}}$, входящем в уравнение каждой моды c индексом $ \pm m$ с определенным коэффициентом усиления ${{g}_{{ \pm m}}}$. Если на исходные управляющие сигналы ${{s}_{n}}(t)$, $n = 1, \ldots ,N$ не накладываются какие-либо ограничения и выбран целевой функционал, который декомпозируется на эти же группы по функциям состояния ${{y}_{{ \pm m}}}$ и управления ${{u}_{k}}$, можно рассматривать $N$ полученных подсистем независимо.

Управляемость этих подсистем связана с темпом убывания коэффициентов ${{g}_{{ \pm m}}}$ и собственных частот ${{\omega }_{{ \pm m}}}$. Доказательство управляемости осложняется еще и тем, что в каждой подсистеме существует две ветви колебаний. При определенном выборе параметров задачи могут совпасть частоты ${{\omega }_{{ - {{m}_{1}}}}}$ и ${{\omega }_{{ + {{m}_{2}}}}}$ для разных номеров ${{m}_{1}} > {{m}_{2}}$ некой группы мод. В этом случае подсистема будет неуправляема из-за существования двух одинаковых осцилляторов при одном управляющем входе.

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

Выберем функцию управления ${{u}_{k}}$ как многочлен от времени четной степени:

где $a_{j}^{k} \in \mathbb{R}$ – неопределенные коэффициенты многочлена. Эти параметры можно использовать, чтобы выполнить, например, ${{N}_{k}}$ терминальных комплексных условий: для нижней ветви в каждой (исключая нулевую) группе мод.

Рассмотрим в качестве примера расчета управляемых движений случай стержня с геометрическими и механическими параметрами (5.1) и тремя пьезоактюаторами ($N = 3$). Поставим задачу гашения первых трех мод колебаний с помощью полиномиального управления (7.1). Для этого нам достаточно выбрать в качестве входных сигналов линейные многочлены:

Терминальные условия (7.2) конкретизируются как

Для демонстрации правильного выбора шести параметров управления (7.3) в достижении предписанного состояния выберем следующее начальное состояние:

Эти условия означают, что первые три моды, которые мы хотим привести в нулевое положение с нулевым импульсом, в момент времени $t = 0$ полностью определяют некоторые упругие перемещения точек системы и нулевую плотность импульса. При этом начальные деформации таковы, что, согласно нормировке (4.12) на каждую моду приходится единичная энергия. Моды с номерами $m = 2Nj$, $j \in \mathbb{N}$, не возбуждаются вовсе. Колебания остальных мод будут определять остаточные движения, энергия которых и покажет качество управления.

Частота самой низшей моды колебаний равна ${{\omega }_{{ - 1}}} = 0.00936$, а соответствующий период ${{T}_{{ - 1}}} = 671$. Чтобы управляющие сигналы не были чрезмерно большими, выберем время управления больше этого периода: $T = 700$. На рис. 6 показаны амплитуды поперечных перемещений центральной линии стержня ${{{v}}_{{ - m}}}(t)$ = $\mathop {\max }\limits_{{{x}_{1}} \in [0,a]} {v}_{2}^{{ - m}}({{x}_{1}},0,0,t)$ первых пяти мод колебаний управляемых, согласно (7.3), линейно. Моды с номерами $m = 1,2,3$, показанные на графике соответственно пунктирной, штрихпунктирной и штриховой линиями точно приходят в нулевое терминальное состояние. Их начальные значения разнятся, так как поперечные отклонения для мод с одинаковой энергией монотонно уменьшаются с увеличением номера. Поведение четвертой и пятой моды, отложенные на рис. 6 длинной штриховой и сплошной линиями соответственно, носят колебательный характер и для выбранного закона управления по энергии сравнимы с приводимыми модами. Для уменьшения остаточных колебаний необходимо повышать степени полиномов управления, используя появившиеся параметры либо для повышения их гладкости на начальном этапе, либо для гашения дополнительных мод.

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

Исследование выполнено за счет гранта Российского фонда фундаментальных исследований (проект № 21-51-12004).