Dynamics of the Euler — Bernoully beam with distributed hysteresis properties

Full Text

Введение

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

Известно, что одним из базовых элементов любой несущей конструкции является балка. Математическое описание балки обычно формализуется в рамках теории Эйлера — Бернулли [40] или теории Тимошенко [41] (рис. 1). Основное различие между вышеозначенными подходами заключается в поведении балки при поперечных деформациях. В теории Эйлера — Бернулли предполагается, что во время деформаций поперечное сечение перпендикулярно нейтральной оси. В то же время в теории Тимошенко оно не перпендикулярно нейтральной оси и напрямую зависит от интенсивности деформации. В настоящей статье модель колебаний балки строится с использованием теории Эйлера — Бернули в предположении, что балка имеет шарнирное закрепление в двух концах и обладает нулевыми начальными условиями (балка покоится и не имеет начальной скорости). Следуя работе [42], классическое уравнение колебаний балки имеет следующий вид:

$\left\{\begin{array}{l}{u}_{tt}+a^2{u}_{xxxx}=g(x,t)t>\mathrm{0,}x\in \left[\mathrm{0,}L\right],\\ u\left(x\mathrm{,0}\right)=u_t\left(x\mathrm{,0}\right)=\mathrm{0,}\\ u\left(\mathrm{0,}t\right)=u(L,t)=S\left(t\right),\\ u_{xx}\left(\mathrm{0,}t\right)=u_{xx}(L,t)=0.\end{array}\right.$ (1)

где $u=u(x,t)$ – отклонение балки от положения равновесия, индексы обозначают частные производные^[1], $g(x,t)$ есть функция, характеризующая вертикальную нагрузку (управление), $S\left(t\right)$ – функция, характеризующая внешнее воздействие, прикладываемое к закрепленным концам балки, $a^2=EI$ – произведение модуля упругости и момента инерции, $t>0$ есть время моделирования, $L$ – длина балки.

Рис. 1. Графическое представление принципиального различия между теориями Тимошенко и Эйлера — Бернулли

Fig. 1. Illustration of the main difference between Timoshenko and Euler — Bernoulli beam theories

Балки в зависимости от своих конструктивных особенностей находят свое применение в различных элементах строительных (и не только) конструкций. Например, двутавровые балки используются для обеспечения устойчивости перекрытий между этажами и в мостовых сооружениях, швеллерные же балки используются для армирования железобетонных конструкций, а также в качестве рам, каркасов и т. д. Кроме того, балки по-разному реагируют на внешнюю нагрузку, которая, в свою очередь, может быть импульсной, вибрационной, подвижной [43], ударной, ветровой. При этом основная задача, решаемая при выборе типа балки применительно к задачам строительства, заключается в обеспечении устойчивости всего сооружения по отношению к потенциальным внешним воздействиям, к числу которых, в первую очередь, относятся сейсмические волны как природного, так и техногенного характера. В таком случае балка выступает в роли демпфирующего элемента, трансформируя входящую энергию воздействия в иные ее формы [44].

Известно, что во время нагрузок балка начинает деформироваться, порождая во внутренней структуре материала силы, стремящиеся вернуть ей исходную форму. Если после снятия нагрузки она принимает форму, отличную от исходной, то следует говорить об упруго-пластическом гистерезисе. Он играет важную роль в современных исследованиях реальных систем и процессов. Это явление было обнаружено во многих областях современной науки: физике [45-53], экономике, биологии, химии и т. д. Недавние исследования показывают важность и целесообразность использования гистерезиса в нейронных сетях, например, включая его как звено или как самостоятельный элемент в функцию активации соответствующей сети. Такой широкий спектр применения обуславливается важными структурными особенностями гистерезиса – зависимостью текущего состояния системы от предыстории и начального состояния как от параметра. Еще более важной особенностью гистерезиса является тот факт, что в ряде задач он позволяет моделировать процесс диссипации энергии, что, в свою очередь, позволяет использовать гистерезисные преобразователи в качестве управляющего элемента, стабилизирующего поведение динамических систем [54-56].

Теоретические модели, описывающие гистерезисные явления, восходят к классическим работам Прандтля [57], Прейзаха [58] и др. В начале 60-х годов XX века появились первые феноменологические модели [59-62]. Теория конструктивных моделей гистерезиса, учитывающая физические особенности конструктивных элементов исследуемых систем, получила свое развитие в трудах М.А. Красносельского и его учеников в Воронежском университете в 60-х годах прошлого века. Фундаментальная монография [63] с подробным описанием свойств конструктивных моделей гистерезиса вышла в 1983 г. Разработанная в ней методология позволила описывать гистерезис на языке нелинейных операторов, зависящих от своего начального состояния как от параметра. Математический аппарат, предложенный в ней, базируется на выделении элементарных носителей гистерезисных свойств – гистеронов – преобразователей, наделенных пространством состояний и входно-выходными соответствиями. Такой подход позволил имплементировать модели гистерезиса в состав более сложных систем, описываемых дифференциальными уравнениями.

В настоящей статье предлагается модель балки с учетом распределенных гистерезисных свойств, в разд. 1 описываются используемые в работе модели гистерезиса, разд. 2 посвящен выводу уравнений колебаний гистерезисной балки, в разд. 3 описаны численные методы, используемые в работе, а в разд. 4 приводятся результаты моделирования и сравнительный анализ классической и гистерезисной балок. В качестве моделируемой нагрузки будут выступать периодические и импульсные воздействия, а также нагрузка сейсмического характера.

1. Модели гистерезиса: феноменологический и конструктивный подходы

1.1. Преобразователь Прандтля

Одной из самых известных и простых гистерезисных моделей является преобразователь Прандтля [64-66] (также известный как “stop operator”). Он используется для моделирования упруго-пластических волокон [67], где состояния полностью определяются величинами переменной деформации и переменного напряжения, а также применяется в качестве элементарной составляющей в сложных континуальных гистерезисных моделях в рамках теории упруго-пластического гистерезиса. Действительно, преобразователь Прандтля идеально описывает поведение материала: материал, находящийся в изначальном состоянии, когда нормальное напряжение и деформация отсутствуют, под нагрузкой начинает вести себя упруго до некоторого порогового значения $h$, а после демонстрирует пластические свойства, находясь под постоянным напряжением. Это поведение также называется “пластическое течение”, так как при уменьшении напряжения материал обладает пластическими свойствами, пока не достигнет порогового значения с противоположным знаком $-h$. Однако отметим, что преобразователь Прандтля позволяет моделировать упруго-пластические свойства без деформационного упрочнения^[2].

Структурно этот преобразователь может быть представлен как последовательное соединение линейного упругого элемента $E$ (например, пружина) с жестким, идеально пластическим элементом $P$ (например, брусок, двигающийся по поверхности, под действием силы сухого трения) (рис. 2). Если предположить, что упругость пружины $E$ равна единице, следует говорить о более простом преобразователе – упоре.

Рис. 2. Графическое представление гистерезисного преобразователя Прандтля. Здесь E и P — упругий и пластический элементы соответственно

Fig. 2. Graphical representation of the Prandtl operator. E and P are the elastic and the plastic elements respectively

Описание преобразователя Прандтля начнем с определения его состояний. Состояниями этого преобразователя являются пары “вход-выход” $\{x,u\}$ с областью определения в виде полосы $\Omega =\Omega ({\Gamma }_P,E)$, угол наклона которой зависит от коэффициента $E$, определяющего модуль упругости. Причем полоса $\Omega$ зависит также и от порогового значения $h$ (полоса пересекает ось $O_x$ в точках $\pm h$), которое является еще одной характеристикой этой модели. Иначе говоря, модель может быть формализована в виде равенства:

$\sigma \left(t\right)={\Gamma }_P\left[h\right]\epsilon \left(t\right),$

 где

${\Gamma }_P\left[h\right]\epsilon \left(t\right)=\left\{\begin{array}{cc}\min \{h,\phi (t,\tau \left)\right\},& если\epsilon \left(t\right)\ge \epsilon (t-\tau ),\\ \max \{-h,\phi (t,\tau \left)\right\},& иначе\end{array}\right.$

$\phi (t,\tau )=E\left[\epsilon \left(t\right)-\epsilon (t-\tau )\right]+\sigma (t-\tau )$

 или в более компактной форме:

$\sigma \left(t\right)={\Gamma }_P\left[h\right]\epsilon \left(t\right)=\min \{h,\max \{-h,\phi (t,\tau )\left\}\right\},$ (2)

 где $\sigma (t-\tau )$ — значение выхода гистерезисного преобразователя в момент времени $t-\tau$, $\epsilon \left(t\right)$ — входной сигнал, $\sigma \left(t\right)$ — выходной сигнал (реакция преобразователя). Графическое представление описанного выше преобразователя отображено на рис. 3.

Рис. 3. График зависимости входа ε(t) от выхода σ(t) для преобразователя Прандтля с пороговым значением, равным h, и областью определения Ω=Ω(Гр, E) в виде полосы с углом наклона, зависящим от коэффициента E

Fig. 3. Dependence of the input ε(t) on the output σ(t) for the Prandtl operator with the threshold value h and the domain Ω=Ω(Гр, E) presented as a band with the slope angle determined by the coefficient E

1.2. Конструктивный подход: преобразователь Прандтля — Ишлинского

Одним из самых известных преобразователей в рамках конструктивного подхода к моделированию упруго-пластического гистерезиса является преобразователь Прандтля — Ишлинского. Обычно этот преобразователь применяется в задачах, где требуется формализовать гистерезисную связь между напряжением и деформацией. Он является компиляцией преобразователей Прандтля и материала Ишлинского. В отличие от преобразователя Прандтля преобразователь Прандтля — Ишлинского позволяет учитывать деформационное упрочнение совместно с упруго-пластическими свойствами (рис. 4). Натурные эксперименты, проведенные с упруго-пластическими материалами, демонстрируют, что состояние материала не только характеризуется входно-выходной характеристикой “напряжение-деформация”, но и зависит от некоторых внутренних переменных. Преобразователь Прандтля — Ишлинского позволяет учитывать предысторию и является удобным инструментом для моделирования гистерезисных зависимостей. Преобразователи с такой структурой относят к классу сложных гистерезисных систем, описываемых посредством континуальных аналогов блок-схем (рис. 5). В настоящем разделе приводится математическое описание этого преобразователя.

Рис. 4. Внутренняя структура континуального преобразователя Прандтля — Ишлинского с пороговыми значениями h=h₁, h₂, ..., h_∞, распределенными с помощью некоторой функции распределения и параллельно соединенными преобразователями Прандтля Гр=Гр₁, Гр₂, ..., Гр_∞

Fig. 4. Inner structure of the Prandtl — Ishlinskii operator with the threshold values h=h₁, h₂, ..., h_∞, distributed using a distribution function and connected in parallel by the Prandtl operators Гр=Гр₁, Гр₂, ..., Гр_∞

Рис. 5. Блок-схема параллельного соединения упоров в модели Прандтля — Ишлинского

Fig. 5. Schema of parallel connection of the stop-operator models in the Prandtl — Ishlinskii model

Зададим семейство однопараметрических преобразователей Прандтля (2), зависящих от порогового значения $h$ как от параметра. Пусть параметр $h$, характеризующий пороговое значение, принадлежит некоторому множеству $\Phi$. Пусть задана неубывающая непрерывная слева функция $\Xi =\Xi \left(h\right)$, $h>0$ на множестве $\Phi$, удовлетворяющая следующим двум условиям:

$\underset{h\to \infty }{}\Xi \left(h\right)=\mathrm{0,}$

${\int }_0^{\infty }\left|\Xi \right(h\left)\right|dh<\infty .$

Также пусть задано множество $Z$ непрерывных функций $\sigma \left(h\right)$, удовлетворяющих условию

$\left|\sigma \right(h\left)\right|\le h.$

Теперь определим множество $\Omega ({\Gamma }_P,E)$ возможных состояний преобразователя ${\Gamma }_P$ как

$\Omega ({\Gamma }_P,E)=\left\{\right\{\epsilon ,{\sigma }_1\left(h\right)\},\{\epsilon ,{\sigma }_2\left(h\right)\left\}\mathrm{,...,}\right\{\epsilon ,{\sigma }_{\infty }\left(h\right)\left\}\right\},$

где $\epsilon$ — это вход, ${\sigma }_i\left(h\right)$ — выход (реакция) $i$-го преобразователя ${{\Gamma }_P}^i$ (каждый преобразователь зависит от порогового значения как от параметра). Входно-выходные соответствия такого преобразователя описываются следующим равенством:

$\sigma \left(t\right)={\Gamma }_{PI}[t_0,{\sigma }_0(h);\Xi ]\epsilon \left(t\right)={\int }_0^{\infty }{\Gamma }_P[t_0,{\sigma }_0(h);h]\epsilon \left(t\right)d\Xi \left(h\right),$ (3)

 а переменное состояние $\left\{\epsilon \right(t),\sigma (h,t\left)\right\}$ определяется равенством

$\sigma (h;t)=Q[t_0,z_0(h);\Xi ]\epsilon \left(t\right)={\Gamma }_P[t_0,{\sigma }_0(h);h]\epsilon \left(t\right),$ (4)

 где $Q$ есть пучок гистеронов (конечный или в общем случае бесконечный набор гистеронов), преобразователь ${\Gamma }_{\text{P}}$ описывает параллельное соединение бесконечного набора преобразователей, описываемых соотношением (2) (см. рис.  4), ${\sigma }_0\left(h\right)$ стоит понимать как значение выхода в предыдущий момент времени (альтернативная запись $\sigma (t-\tau ,h)$).

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

$\sigma \left(t\right)=\sum _{i=1}^{N}f\left(h_i\right){\Gamma }_P[t_0,{\sigma }_i;h_i]\epsilon \left(t\right).$ (5)

Большое разнообразие работ, использующих преобразователь Прандтля — Ишлинского, можно условно разделить на две категории: фундаментальные и прикладные. В рамках фундаментальных работ исследуются структурные особенности преобразователя [68-70], вопросы корректности, а также проблемы идентификации его параметров [71]. В работах прикладной направленности преобразователь используется для описания гистерезисных явлений в механических системах, например актуаторы, магнитореологические эластомеры [72], несущие конструкции [73] и т. д.

1.3. Феноменологический подход: преобразователь Боука — Вена

Наряду с конструктивным, не менее известным является феноменологический подход к моделированию гистерезиса. Он позволяет моделировать петли гистерезиса разной структуры, не принимая во внимание физическую основу моделируемого процесса. Обычно феноменологический подход применяется в том случае, когда необходимо оценить реакцию системы на сигнал лишь на основе экспериментальных данных о характеристиках входного воздействия и отклике системы на него. Одним из самых известных преобразователей в рамках этого подхода является преобразователь Боука — Вена, предложенный в [21; 22]. Входные-выходные соответствия преобразователя Боука — Вена формализуются уравнением:

${\dot{\Gamma }}_{BW}=\dot{x}\left\{A-\left[\beta sign\left(\dot{x}{\Gamma }_{BW}\right)+\gamma \right]|{\Gamma }_{BW}{|}^n\right\},$ (6)

где $x=x\left(t\right)$ — входной сигнал, ${\Gamma }_{BW}={\Gamma }_{BW}\left(t\right)$ — выходной сигнал (реакция преобразователя), $\gamma$, $n$, $\beta$, $A>0$ — безразмерные параметры. Как можно заметить, преобразователь Боука — Вена является четырехпараметрическим, что позволяет моделировать различные по форме и площади гистерезисные петли. Отметим, что проблеме идентификации этих параметров посвящено достаточно большое количество публикаций.

Очевидным преимуществом преобразователя Боука — Вена являются простота численной реализации и развитый математический аппарат [74-77]. В этих работах рассмотрены задачи возникновения предельного цикла, взаимная зависимость параметров преобразователя, адаптация модели применительно к некоторому набору прикладных задач и т. д. На сегодняшний день преобразователь Боука — Вена применяется и для моделирования процессов различного масштаба [15; 39].

2. Уравнение колебаний гистерезисной балки:
гамильтоновский подход

 В настоящем разделе производится вывод уравнений колебаний балки с распределенными гистерезисными свойствами (рис. 6). Отметим, что в подавляющем большинстве прикладных задач рассматривались объекты с сосредоточенными гистерезисными свойствами. Распределенные гистерезисные свойства рассматривались в весьма ограниченном количестве работ, из которых отметим работу [73]. Ниже приводится вывод уравнений колебаний гистерезисной балки, основанный на вариационном принципе, а именно принципе наименьшего действия Гамильтона. В рамках этого подхода искомые уравнения являются следствием экстремального значения функционала действия.

Рис. 6. Изгиб поперечного сечения балки во время деформации

Fig. 6. Bending of the cross section of the beam during the deformation

Рассмотрим стержень, работающий на изгиб, имеющий постоянную площадь поперечного сечения $S=const$, симметричную относительно обеих осей $O_y$ и $O_z$, при этом нейтральная ось стержня совпадает с осью $O_x$. Отметим, что в текущем выводе учитываются только поперечные колебания в плоскости $O_{xy}$, продольными колебаниями пренебрегают.

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

$I\left(u\right)={\int }_{t_0}^{t_1}Ldt,$ (7)

где $t_0\le t\le t_1$ — некоторый промежуток времени, $u=u(x,t)$ — кривая, описывающая динамику балки, $L$ — функция Лагранжа, равная разности кинетической и потенциальной энергии:

$L=T-U,$ (8)

где $U$ — потенциальная энергия, $T$ — кинетическая энергия. Определим потенциальную энергию балки как

$U=\frac{1}{2}{\int }_0^{l}M\chi dx,$ (9)

 где $M$ — изгибающий момент, а $\chi$ — деформация балки в точке $x$. Деформацию определим через кривизну.

Пусть кривая, описывающая положение балки в плоскости $O_{xy}$, задана уравнением $F=F(x,y)=0$ ($F\in {ℝ}^2$), тогда ее кривизна будет описываться уравнением:

$\kappa =\frac{|F_y^2{F}_{xx}-F_x{F}_y{F}_{xy}+F_x^2{F}_{yy}|}{{(F_x^2+F_y^2)}^{3/2}}.$

 Так как продольными колебаниями можно пренебречь (согласно предположению), то уравнение, описывающее положение балки, будет иметь вид $F(x,y)=y-u\left(x\right)=0$, а ее кривизна определяться следующим соотношением:

$\kappa =\frac{u_{xx}}{{(1-u_x^2)}^{3/2}}.$

Ограничиваясь рассмотрением лишь малых углов отклонения поперечного сечения $\varphi \ll 1$ ($u_x\ll 1$ и как следствие $1-u_x^2\underset{u_x\to 0}{\to }1$), кривизна будет определяться равенством:

$\kappa =u_{xx}.$

 В свою очередь, изменение кривизны балки равно^[3]

$\chi =k-k_0=k-0=k=u_{xx}$ (10)

 где $k_0$ — начальная кривизна покоящейся балки, равная нулю.

Для вывода уравнения изгибающего момента выберем два параллельных сечения $\Omega x$ и $\Omega (x+\Delta x)$ на расстоянии $∆x$. Очевидно, что углы поворота в этих сечениях будут различными во время деформации. Вследствие этого происходят растяжение (в верхней части балки) и сжатие (в нижней части балки) материала балки в продольном направлении. Рассмотрим часть балки, находящейся на расстоянии $r=r\left(y\right)$ от ее оси симметрии. Как видно на рис. 7, длина растянутого материала изменилась на $\Delta s=r\Delta \varphi$. Относительное удлинение материала

$\epsilon (x,\Delta x,y,t)=-\frac{\Delta s}{\Delta x}=-r\frac{\Delta \varphi }{\Delta x},$ (11)

 а деформация растяжения-сжатия в точке $x$

$\epsilon (x,y,t)=-\underset{\Delta x\to 0}{}\epsilon (x,\Delta x,y,t)=-r\underset{\Delta x\to 0}{}\frac{\Delta \varphi }{\Delta x}=-r{\varphi }_x.$ (12)

Рис. 7. Деформированная часть балки, полученная двумя параллельными сечениями Ω(x) и Ω(x+Δx)

Fig. 7. Deformed part of the beam obtained by means of two parallel cross sections Ω(x) and Ω(x+Δx)

 Угол отклонения поперечного сечения  от “нормального” положения в точке $x$ (см. рис. 6) с учетом ортогональности оси симметрии направлению деформации равен:

$\varphi =u_x\mathrm{.}$ (13)

 Таким образом, подставляя (13) в (14), получим

$\epsilon (x,y,t)=-r{u}_{xx}.$ (14)

Изгибающий момент, возникающий относительно оси $O_z$ вследствие нормального напряжения в сечении $\Omega \left(x\right)$, описывается уравнением

$M(x,t)={\int }_{\Omega }r\sigma dS.$ (15)

Теперь вычислим нормальное напряжение в сечении $\Omega \left(x\right)$, полагая, что материал балки обладает упруго-пластическими свойствами с деформационным упрочнением. Для этого воспользуемся преобразователем Прандтля — Ишлинского:

$\sigma \left(t\right)={\Gamma }_{PI}[t_0,{\sigma }_0;\Xi ]\epsilon \left(t\right):={\int }_0^{\infty }{\Gamma }_P[t_0,{\sigma }_0;h]\epsilon \left(t\right)d\Xi \left(h\right),$ (16)

где $\sigma \left(t\right)$ — выход (реакция преобразователя), ${\Gamma }_{PI}$ — преобразователь Прандтля — Ишлинского, $\epsilon \left(t\right)$ — вход. Выполним переход от одномерного преобразователя Прандтля — Ишлинского к двумерному, введя пространственный параметр $\zeta$, принадлежащий множеству измеримых параметров, по Лебегу, как предлагается в работе [73]. Тогда параметрический преобразователь примет вид:

$\sigma (\zeta ,t)={\Gamma }_{PI}[t_0,{\sigma }_0;\Xi ]\epsilon (\zeta ,t):={\int }_0^{\infty }{\Gamma }_P[t_0,{\sigma }_0;h]\epsilon (\zeta ,t)d\Xi (\zeta ,h),$

где $\zeta$ — параметр, характеризующий пространственную (по материалу балки) плотность распределения элементарных носителей гистерезиса. Или, используя следующее соотношение $\zeta =(x,y)$ (параметр $\zeta$ представлен в виде пространственной переменной), окончательно получим:

$\sigma (x,y,t)={\Gamma }_{PI}[t_0,{\sigma }_0;\Xi ]\epsilon (x,y,t):={\int }_0^{\infty }{\Gamma }_P[t_0,{\sigma }_0;h]\epsilon (x,y,t)d\Xi (x,y,h).$ (17)

 С учетом соотношения (17) уравнение (15) примет вид:

$M(x,t)={\int }_{\Omega }r\left({\int }_0^{\infty }{\Gamma }_P[t_0,{\sigma }_0;h]\epsilon (x,y,t)d\Xi (x,y,h)\right)dS.$ (18)

 Раскрыв скобки и учитывая соотношение (14), получим

$M(x,t)={\int }_{\Omega }{\int }_0^{\infty }r{\Gamma }_P[t_0,{\sigma }_0;h](-r{u}_{xx})d\Xi (x,y,h)dS.$ (19)

Теперь воспользуемся свойством преобразователя Прандтля [73], а именно соотношением

$r{\Gamma }_P\left[h\right](-r{u}_{xx})=-r^2\Gamma [h/|r\left|\right]u_{xx}.$ (20)

Для подынтегрального выражения из соотношения (19) с учетом соотношения (14)

$r{\Gamma }_P\left[h\right](-r{u}_{xx})=-r^2\Gamma [h/|r\left|\right]u_{xx}.$ (21)

и введя переменную $\widehat{\Xi }(x,h)={\int }_{\Omega }{r}^2\Xi (x,y,|r\left|h\right)dS$, получим

$M(x,t)=-{\int }_0^{\infty }{\Gamma }_P[t_0,{\sigma }_0;h]u_{xx}{d}_h\widehat{\Xi }(x,h).$ (22)

Как можно заметить, уравнение (22) совпадает с уравнением (16), что позволяет записать соотношение в следующей форме:

$M(x,t)=-{\Gamma }_{PI}{u}_{xx}.$ (23)

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

Окончательно принимая во внимание соотношения (10) и (23), уравнение для потенциальной энергии примет вид

$U=-\frac{1}{2}{\int }_0^l{u}_{xx}{\Gamma }_{PI}\left[u_{xx}\right]dx.$ (24)

 Кинетическую энергию определим классическим образом:

$T=\frac{1}{2}{\int }_0^l\rho \left(x\right)S{u}_t^{2}dx.$ (25)

 Тогда интеграл действия для такой системы примет вид

$I\left(u\right)={\int }_{t_1}^{t_2}\frac{1}{2}\left({\int }_0^L{u}_{xx}{\Gamma }_{PI}\left[u_{xx}\right]dx+{\int }_0^L\rho S{u}_t^{2}dx\right)dt.$ (26)

 Из принципа наименьшего действия в форме Гамильтона следует, что функция $u=u(x,t)$ является экстремалью функционала (26). Рассмотрим вариацию этой функции $\widehat{u}(x,t)=u(x,t)+\epsilon \eta (x,t)$ , для которой справедливо следующее: 

• в начальный и конечный моменты времени функция $u\left(x,t\right)$ не варьируется;
• $\eta (x,t_1)=\eta (x,t_2)=0$ для таких $x$, для которых верно $0\le x\le L$.

 Таким образом, функционал (26) для кривой $\hat{u}$ примет вид

$I={\int }_{t_1}^{t_2}\frac{1}{2}\left({\int }_0^L(u_{xx}+\int {\eta }_{xx})\Gamma [u_{xx}+\int {\eta }_{xx}]dx+\rho S{\int }_0^L{(u_t+\int {\eta }_t)}^{2}dx\right)dt.$ (27)

 Ясно, что функционал действия (27) достигает экстремума в точке $\epsilon =0$. Согласно необходимому условию существования экстремума, положим вариацию функционала действия относительно $\epsilon$ равной нулю

${\left.\frac{dI}{d\int }\right|}_{\int =0}=0.$ (28)

Перед вычислением производной по $\epsilon$ раскроем скобки внутри подынтегральных уравнений (27):

${\int }_0^L{u}_{xx}{\Gamma }_{PI}[u_{xx}+\epsilon {\eta }_{xx}]+\epsilon {\eta }_{xx}{\Gamma }_{PI}[u_{xx}+\epsilon {\eta }_{xx}]dx,$ (29)

$\rho S{\int }_0^L(u_t^2+2\epsilon u_t{\eta }_t+{\eta }_t^2)dx.$ (30)

 Подставляя в (28) выражение интеграла действия (27) с учетом (29) и (30), вычислим производную по $\epsilon$:

${\int }_{t_1}^{t_2}{\left((I_1+I_2)+\rho S{I}_3\right)}_{\epsilon }dt=\mathrm{0,}$ (31)

 где

$I_1={\int }_0^L{u}_{xx}{\Gamma }_{PI}[u_{xx}+\epsilon {\eta }_{xx}]dx,$ (32)

$I_2={\int }_0^L\epsilon {\eta }_{xx}{\Gamma }_{PI}[u_{xx}+\epsilon {\eta }_{xx}]dx,$ (33)

$I_3={\int }_0^L(u_t^2+2\epsilon u_t{\eta }_t+{\eta }_t^2)dx.$ (34)

Предположим, что функция $\eta (x,t)$ на концах равняется нулю, то есть $\eta \left(\mathrm{0,}t\right)=\eta (L,t)=0$. Рассмотрим более детально каждое слагаемое (32)–(34).

Найдем частную производную по $\epsilon$ для (32). Для этого продифференцируем подынтегральное слагаемое^[4]:

${\int }_0^L{u}_{xx}{\left({\Gamma }_{PI}[u_{xx}+\epsilon {\eta }_{xx}]\right)}_{\epsilon }dx.$ (35)

 Для дальнейших вычислений предположим, что весовая функция $\widehat{\Xi }$ дифференцируема и ее производная есть функция $f\left(h\right)$. Следуя этому, преобразователь Прандтля — Ишлинского примет вид:

${\Gamma }_{PI}\left[h\right]\int ={\int }_0^{\infty }f\left(h\right){\Gamma }_P\left[h\right]\int dh.$ (36)

 Учитывая (36), запишем полную форму производной преобразователя Прандтля — Ишлинского из соотношения (35):

${\left({\Gamma }_{PI}[u_{xx}+\epsilon {\eta }_{xx}]\right)}_{\epsilon }={\int }_0^{\infty }f\left(h\right){\left({\Gamma }_P[u_{xx}+\epsilon {\eta }_{xx}]\right)}_{\epsilon }dh,$ (37)

где ${\Gamma }_{\text{P}}$ — преобразователь Прандтля. Представляя преобразователь Прандтля в форме (2):

${\Gamma }_P\left[\varphi \right(x,t,\epsilon \left)\right]=\min \{h,\max \{-h,\psi (x,t,t_0,\epsilon )\left\}\right\},$ (38)

 где $\psi (x,t,t_0,\epsilon )={\Gamma }_{P0}+\varphi (x,t,\epsilon )-\varphi (x,t_0,\epsilon )$ и $\varphi (x,t,\epsilon )=u_{xx}+\epsilon {\eta }_{xx}$, его производная примет вид:

${\Gamma }_P{\left[\varphi \right(x,t,\epsilon \left)\right]}_{\epsilon }=\left\{\begin{array}{l}\mathrm{0,}если\psi (x,t,t_0,\epsilon )\in [-h,h],\\ {\eta }_{xx}(x,t)-{\eta }_{xx}(x,t_0),иначе.\end{array}\right.$ (39)

 Полученный выше результат можно представить в виде

${\Gamma }_P{\left[\varphi \right(x,t,\epsilon \left)\right]}_{\epsilon }=\left({\eta }_{xx}\right(x,t)-{\eta }_{xx}(x,t_0\left)\right)\times$

$\times \left(\Theta \right(\psi (x,t,t_0,\epsilon )+h)-\Theta (\psi (x,t,t_0,\epsilon )-h\left)\right),$ (40)

 где $\Theta (\cdot )$ — функция Хевисайда. Тогда, подставляя полученный результат в (37) и вынося слагаемые, не зависящие от $h$, получим

${\left({\Gamma }_{PI}[u_{xx}+\epsilon {\eta }_{xx}]\right)}_{\epsilon }=\left({\eta }_{xx}\right(x,t)-{\eta }_{xx}(x,t_0\left)\right)\times$

$\times {\int }_0^{\infty }f\left(h\right)\left(\Theta \right(\psi (x,t,t_0,\epsilon )+h)-\Theta (\psi (x,t,t_0,\epsilon )-h\left)\right)dh.$ (41)

 Подставив в (35) выражение (41), обозначив интеграл как некоторую функцию $\beta (x,t,t_0,\epsilon )$:

${\int }_0^L{u}_{xx}\left({\eta }_{xx}\right(x,t)-{\eta }_{xx}(x,t_0\left)\right)\beta (x,t,t_0,\epsilon )dx,$ (42)

 проинтегрируем два раза по частям, принимая за $v$ функцию от $\eta$ (чтобы избавиться от двух производных по $x$), и окончательно получим

${\left(I_1\right)}_{\epsilon }={\int }_0^L\left(\eta \right(x,t)-\eta (x,t_0\left)\right)\alpha (x,t,t_0,\epsilon )dx,$ (43)

 где $\alpha (x,t,t_0,\epsilon )$ — некоторая произвольная функция, полученная в ходе интегрирования.

Теперь рассмотрим второй интеграл $I_2$. Дифференцируя интеграл по $\epsilon$:

${\int }_0^L{\eta }_{xx}{\Gamma }_{PI}[u_{xx}+\epsilon {\eta }_{xx}]dx+{\int }_0^L\epsilon {\eta }_{xx}{\left({\Gamma }_{PI}\right[u_{xx}+\epsilon {\eta }_{xx}\left]\right)}_{\epsilon }dx$ (44)

 и полагая $\epsilon =0$, получим следующее выражение:

${\int }_0^L{\eta }_{xx}{\Gamma }_{PI}\left[u_{xx}\right]dx.$ (45)

 Теперь продифференцируем выражение два раза по частям и окончательно получим

${\left(I_2\right)}_{\epsilon }={\int }_0^L\eta {\Gamma }_{PI}{\left[u_{xx}\right]}_{xx}dx.$ (46)

Нетрудно заметить, что в интеграле (34) только одно слагаемое зависит от $\epsilon$. Следовательно, дифференцируя по $\epsilon$, имеем:

${\int }_0^{L}2u_t{\eta }_{t}dx.$ (47)

 Теперь, продифференцировав один раз по частям полученное выражение, окончательно имеем:

${\left(I_3\right)}_{\epsilon }={\int }_0^{L}2\eta u_{tt}dx.$ (48)

Таким образом, подставляя в (31) выражения (43), (46), (48), получим:

${\int }_{t_1}^{t_2}\left({\int }_0^L\left(\eta \right(x,t)-\eta (x,t_0\left)\right)\alpha (x,t,t_0\mathrm{,0})dx+{\int }_0^L\eta \Gamma {\left[u_{xx}\right]}_{xx}dx\right)+$

$+\rho S{u}_{tt}\eta dxdt=0.$

(49)

 Заметим, что вышеполученное уравнение есть необходимое условие Эйлера локального экстремума функционала действия. Воспользуемся леммой Дебуа-Реймона, которая является основой для построения теории обобщенных функций, применительно к уравнению (49) и получим окончательно уравнение колебаний гистерезисной балки^[5]:

$\Gamma {\left[u_{xx}\right]}_{xx}+\rho S{u}_{tt}=0.$ (50)

Учитывая потенциал внешней нагрузки и полагая произведение распределением массы на площадь поперечного сечения, равного единице, перепишем уравнение в виде

$u_{tt}+\Gamma {\left[u_{xx}\right]}_{xx}=g(x,t),$ (51)

 где $g(x,t)$ — функция внешней нагрузки.

3. Численные методы: явная разностная схема и метод Рунге — Кутты

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

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

Введем в области $\Omega =\{0\le x\le L\mathrm{,0}\le t\le T\}$ равномерную сетку с шагом по пространству $h_x$ и с шагом по времени $h_t$, равную

$h_x=\frac{L}{N}{h}_t=\frac{T}{M},$ (52)

 где $M$ — количество узлов сетки по временной координате, $T$ — время моделирования, $N$ — количество узлов сетки по пространственной координате (рис. 8). Причем шаги подчиняются соотношению $h_t\ll h_x$ для устойчивости описываемой схемы. Первые два временных шага $t=0$, $t=1$ вычисляются на основе начальных условий. Используя граничные условия, нетрудно идентифицировать все значения в узлах с координатами $(x_i,t_j)$, где $i=\mathrm{0,1,}N-\mathrm{1,}N-2$, $j=\overline{\mathrm{0,}M}$. Для вычисления остальных значений аппроксимируем уравнение колебаний гистерезисной балки.

Рис. 8. Область $\Omega =\{0\le x\le L\mathrm{,0}\le t\le T\}$ с введенной равномерной сеткой размера M на N. Узлы сетки, обозначенные белым цветом, вычисляются в ходе решения уравнения, узлы, обозначенные черным цветом, вычисляются из граничных условий

Fig. 8. Grid $\Omega =\{0\le x\le L\mathrm{,0}\le t\le T\}$ with size M for N. Undefiined grid nodes are denoted by white color, nodes denoted by black color are calculating based on boundary conditions

Уравнения, аппроксимирующие вторые производные по времени и пространству, примут вид:

$u_{tt}=\frac{u_i^{j-1}-2u_i^j+u_i^{j+1}}{h_t^2},$ (53)

$u_{xx}=\frac{u_{i-1}^j-2u_i^j+u_{i+1}^j}{h_x^2},$ (54)

 где $u_i^j$ соответствует значению в области $\Omega$ с координатами $\left(i,j\right)$. Аппроксимация нулевых начальных и граничных условий очевидна, поэтому не приводится в настоящей статье. Аппроксимация функции нагрузки $g(x,t)$ выполнена следующим образом:

$g(x,t)=g(j{h}_x,i{h}_t),$ (55)

 а аппроксимация гистерезисного оператора имеет вид:

$\Gamma {\left[u_{xx}\right]}_{xx}=\frac{\Gamma {\left[u_{xx}\right]}_i^{j-1}-2\Gamma {\left[u_{xx}\right]}_i^j+\Gamma {\left[u_{xx}\right]}_i^{j+1}}{h_x^2},$ (56)

 и, вводя новую переменную ${\tilde{u}}_i^j=\frac{u_{i-1}^j-2u_i^j+u_{i+1}^j}{h_x^2}$, согласно (54), окночательно получим соотношение:

$\Gamma \left[u_{xx}\right]=\frac{\Gamma {\left[{\tilde{u}}_i^j\right]}_i^{j-1}-2\Gamma {\left[{\tilde{u}}_i^j\right]}_i^j+\Gamma {\left[{\tilde{u}}_i^j\right]}_i^{j+1}}{h_x^2},$ (57)

где $\Gamma$ — гистерезисный оператор. Подставляя полученные аппроксимации в уравнение (51), получим следующее соотношение:

$\frac{u_i^{j-1}-2u_i^j+u_i^{j+1}}{h_t^2}+\frac{\Gamma {\left[{\tilde{u}}_i^j\right]}_i^{j-1}-2\Gamma {\left[{\tilde{u}}_i^j\right]}_i^j+\Gamma {\left[{\tilde{u}}_i^j\right]}_i^{j+1}}{h_x^2}=g(j{h}_x,i{h}_t).$ (58)

 Оставляя в левой части уравнения только неизвестные, окончательно получим уравнение колебаний гистерезисной балки, аппроксимированное явной численной схемой:

$u_i^{j+1}=2u_i^j+-u_i^{j-1}{h}_t^{2}g(j{h}_x,i{h}_t)-h_t^2\frac{\Gamma {\left[{\tilde{u}}_i^j\right]}_i^{j-1}-2\Gamma {\left[{\tilde{u}}_i^j\right]}_i^j+\Gamma {\left[{\tilde{u}}_i^j\right]}_i^{j+1}}{h_x^2}.$ (59)

 Таким образом, значение на $j+1$ временном слое будет высчитываться на основе $j$ и $j-1$ слоев.

В зависимости от формализации гистерезисного звена будут рассматриваться различные подходы к его разностностной дискретизации. Так, для расчета реакции преобразователя Прандтля — Ишлинского будет использоваться его дискретный аналог (5). Выход преобразователя Боука — Вена рассчитывается, используя модификацию хорошо известного метода Рунге — Кутты четвертого порядка. Кратко опишем используемый алгоритм.

Пусть $F=F(\xi ,\eta )$ есть правая часть дифференциального уравнения преобразователя Боука — Вена со входом $\xi$ и значением этой функции в предыдущий момент времени $\eta$. Введем вектор

$F^j=\left(\begin{array}{c}F({\xi }_2^j,{\eta }_2^j)\\ F({\xi }_3^j,{\eta }_3^j)\\ ⋮\\ F({\xi }_{N-2}^j,{\eta }_{N-2}^j)\end{array}\right)$ (60)

 преобразователей Боука — Вена со значениями ${\xi }_i^j=\frac{u_{i-1}^j-2u_i^j+u_{i+1}^j}{h_x^2}$ и ${\eta }_i^j={\Gamma }_{i+1}^{j+1}$ соответственно. Тогда модель Рунге — Кутты 4-го порядка применительно к решению уравнения Боука — Вена можно записать в векторном виде:

${\Gamma }^{j+1}={\Gamma }^j+\frac{1}{6}({k_1}^j+2{k_2}^j+2{k_3}^j+{k_4}^j),$ (61)

где

${k_1}^j=h_t{F}^j{|}_{{\xi }_i^j={\xi }_i^j},$

${k_2}^j=h_t{F}^j{|}_{{\xi }_i^j={\xi }_i^j+{k_1}_i^j/2},$

${k_3}^j=h_t{F}^j{|}_{{\xi }_i^j={\xi }_i^j+{k_2}_i^j/2},$

${k_4}^j=h_t{F}^j{|}_{{\xi }_i^j={\xi }_i^j+{k_3}_i^j},$

 где $F^j$, ${k_1}^j, {k_2}^j, {k_3}^j, {k_4}^j,$есть векторы-столбцы на  временном интервале. Таким образом, представленный метод позволяет оптимизировать вычисления, значительно сократив время их выполнения (в средах, где есть поддержка векторного исчисления, например, среда Mathlab или библиотека NumPy для языка Python и др.).

4. Результаты моделирования колебаний

 В настоящем разделе описываются результаты моделирования колебаний классической и гистерезисной балок в зависимости от разных типов нагрузок. Моделирование феномена гистерезиса осуществляется посредством двух преобразователей: Боука — Вена и Прандтля — Ишлинского. Также в разделе приводится сравнительная характеристика между этими моделями. Начальные условия полагаются нулевыми (балка покоится в начальный момент времени и не имеет начальной скорости соответственно) и ненулевыми граничные условия, когда моделируется сейсмическая нагрузка (на концах моделируется воздействие сейсмической волны) и нулевыми для иных нагрузок (рис. 9). Характеристики балки постоянны во всех экспериментах (кроме тех случаев, где явно указаны другие параметры): 

• длина балки балки равна $L=1$;
• время моделирования равно $T=5$;
• коэффициент $a^2$ равняется 0.5;

 параметры численного моделирования полагались следующими: 

• шаг по времени равен $h_t=0.00125$;
• шаг по пространсву равен $h_x=0.05$;

 а параметры гистерезисных моделей: 

• пороговые значения $h_i$ преобразователя Прандтля — Ишлинского распределялись в интервале $\left(0,100\right)$ с шагом $0,5$;
• коэффициент $E$ преобразователя Прандтля — Ишлинского полагается равным 0.9.
• параметры преобразователя Боука — Вена полагаются равными: $A=0.9, \beta =0.2, \gamma =0.5, n=0.2$.

Рис. 9. Шарнирное закрепление балки с учетом внешней нагрузки g(x,t) (управления) и сейсмологической нагрузки S(t) в месте закрепления

Fig. 9. Fixed boundary of the beam taking into account the external load g(x,t) (control) and seismological load S(t) at the boundary connection  

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

4.1. Моделирование нагрузки как функции периодической структуры

Одним из классических типов нагрузки является периодическая нагрузка. С ее помощью можно понять, как реагирует система на периодические воздействия с меняющейся частотой или амплитудой. Обычно в качестве периодической нагрузки принято использовать тригонометрические функции, изменяющиеся в зависимости от времени, вида $A\sin (\omega t+{\varphi }_0)$ или $A\cos (\omega t+{\varphi }_0)$, где $A$ — амплитуда, $\omega$ — частота и ${\varphi }_0$ — начальная фаза колебаний. В настоящей статье в качестве функции периодического воздействия будет использоваться функция $g(x,t)=0.12\sin (2\pi /2-4\pi t)$ (рис. 10).

Рис. 10. График функции нагрузки $g(x,t)=0.12\sin (2\pi /2-4\pi t)$

Fig. 10. Graphical representation of the load function $g(x,t)=0.12\sin (2\pi /2-4\pi t)$

Как можно видеть из полученных результатов (рис. 11), реакция каждого типа балки на периодическое воздействие разная. Классическая балка с начального момента времени начинает демонстрировать квазипериодические формы колебаний c двумя частотами (основной с наибольшей мощностью, и второстепенная, с меньшей), что подтверждается анализом спектра Фурье-преобразования (рис. 12, а). Балка с гистерезисными свойствами, формализуемыми преобразователем Боука — Вена, демонстрирует схожие результаты с классической балкой. А именно на начальном временном отрезке можно видеть схожие колебания с классической моделью. Однако, начиная с некоторого момента времени, устанавливаются периодические колебания с постоянной частотой (отметим, что частота у каждого из преобразователей отличается, ровно как и амплитуда колебаний). Анализ спектра Фурье-преобразования показывает, что использование преобразователя Боука — Вена способствует уменьшению мощности неосновной частоты колебаний (рис. 12, в) практически нивелируя ее. Используя в качестве носителя гистерезисных свойств преобразователь Прандтля — Ишлинского, можно добиться лучших результатов в сравнении с остальными моделями. А именно с ее помощью, во-первых, можно добиться уменьшения амплитуды основной частоты, во-вторых, редуцировать второстепенную частоту (рис. 12, б). Представленные гистерезисные зависимости (рис. 11, в, г) позволяют оценить энергетический вклад гистерезисного звена в редуцирование колебаний: площадь петли пропорциональна энергии, “отводимой” гистерезисным звеном из системы.

Рис.11. Графики колебаний балки в точке x=L/2 в зависимости от времени: синяя кривая — классическая модель балки, оранжевая кривая — гистерезисная балка с использованием преобразователя Прандтля — Ишлинского, красная кривая — гистерезисная балка с использованием преобразователя Боука — Вена ( а, б); параметрические графики гистерезисных петель (время t как параметр): входно-выходная зависимость входа u_xx от выхода Γ[u_xx] в точке x=L/2 ( в, г)

Fig. 11.Vibration graphs at the point x=L/2 depending on time: the blue curve denotes the classical beam, the orange curve denotes the Prandtl — Ishlinskii model, the red curve denotes the Bouc — Wen model ( а, b); parametric graphs of hysteresis loops (time t as a parameter): input-output dependence of the input u_xx on the output Γ[u_xx] at the point x=L/2 ( c, d) 

Рис. 12. Спектр Фурье-преобразования колебаний балки в точке x=L/2: а — классическая модель, б — гистерезисная балка с использованием модели Прандтля — Ишлинского, в — гистерезисная балка с использованием модели Боука — Вена

Fig. 12. Fourier spectrum of the vibrations of the beam at the point x=L/2: left — classical model, middle — Prandtl — Ishlinskii model, right — Bouc — Wen model 

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

4.2. Моделирование воздействия импульсной нагрузки

 Нередко в рамках прикладных задач важно знать реакцию несущей конструкции на импульсную нагрузку. Самым простым примером такой нагрузки является точечный удар по конструкции с некоторой “силой”. Классическим методом формализации импульсной нагрузки является дельта-функция Дирака. Однако, принимая во внимание использование численных методов в работе, имеет смысл рассмотреть подобную ей функцию, которая в одномерном случае имеет вид:

$\widehat{\delta }\left(t\right)=\left\{\begin{array}{cc}A,& еслиt=\mathrm{0,}\\ \mathrm{0,}& иначе,\end{array}\right.$

 а в двумерном:

$\widehat{\delta }(x,t)=\left\{\begin{array}{cc}A,& еслиt=0иx=\mathrm{0,}\\ \mathrm{0,}& иначе,\end{array}\right.$

 где $A$ — “сила” удара. Использование такой функции позволяет нивелировать бесконечности в расчетах, что сильно упрощает задачу моделирования. Ясно, что балка за счет внутренних сил должна “поглощать” удар, преобразуя энергию удара в иные ее формы, например, в тепловую энергею (нагрев балки). В настоящей статье исследуется реакция классической и “гистерезисной” балок на импульсную нагрузку (удар), приложенную к ее определенной дискретной точке в некоторый момент времени. Таким образом, функцию $g(x,t)$ можно представить в виде:

$g(x,t)=\widehat{\delta }(x-\chi ,t-\tau )=\widehat{\delta }(x-\chi )\widehat{\delta }(t-\tau ),$

 где $\chi \in \left[\mathrm{0,}L\right]$ есть дискретная точка балки, в которой действует импульсная нагрузка, а $\tau \in \left(\mathrm{0,}T\right)$ — момент времени, в который действует нагрузка.

Результаты моделирования, продемонстрированные на рис. 13, показали, что учет гистерезисных свойств приводит к совершенно иной динамике колебаний, нежели в классической модели. А именно начиная с момента времени $t=\tau$ классическая балка начинает совершать колебания с постоянной частотой, причем в ходе экспериментов было установлено, что импульсная нагрузка любой “силы” выводит балку из состояния равновесия, заставляя балку бесконечно долго колебаться. Напротив, балка с учетом гистерезисных свойств с момента времени  начинает отклоняться от положения равновесия, достигает некоторой точки, совершая при этом одиночное колебание малой амплитуды, после чего достигает устойчивого положения, отличного от положения равновесия и остается в этом положении бесконечно долго. При этом на рис. 13,  в можно видеть, что результаты, демонстрируемые с использованием преобразователей Прандтля — Ишлинского и Боука — Вена, отличаются, однако схожи по своей структуре. Схожесть результатов подтверждается и подобностью гистерезисных петель, продемонстрированных на рис. 14.

Рис. 13. Графики колебаний балки в точке x=L/2 в зависимости от времени при импульсной нагрузке с силой A=1000 и параметрами χ=0.5, τ=0.5: синяя кривая — классическая модель балки, оранжевая кривая — гистерезисная балка с использованием преобразователя Прандтля — Ишлинского, красная кривая — гистерезисная балка с использованием преобразователя Боука — Вена

Fig. 13. Vibration graphs at the point x=L/2 depending on time under the impulse load with the force A=1000 and parameters χ=0.5, τ=0.5: the blue curve denotes the classical beam, the orange curve denotes the Prandtl — Ishlinskii model, the red curve denotes the Bouc — Wen model 

Рис. 14. Параметрические графики гистерезисных петель (время t как параметр): входно-выходная зависимость входа u_xx от выхода Γ[u_xx] в точке x=L/2: а — преобразователь Прандтля — Ишлинского, б — преобразователь Боука — Вена

Fig. 14. Parametric graphs of hysteresis loops (time t as a parameter): input-output dependence of the input u_xx on the output Γ[u_xx] at the point x=L/2: a — the Prandtl — Ishlinskii operator, b — the Bouc — Wen operator 

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

Рис. 15. Красной кривой соответствует положение гистерезисной балки, формализуемой преобразователем Боука — Вена, синей кривой — положение классической балки, зеленая линия демонстрирует место, в которое был совершен удар

Fig. 15. The red curve corresponds to the position of the hysteretic beam formalised using the Bouc — Wen operator, the blue curve corresponds to the position of the classical beam, and the green curve corresponds to the point of impact  

4.3. Моделирование воздействия сейсмической нагрузки

 В настоящем разделе приводятся результаты моделирования колебаний классической и “гистерезисной” балок, находящихся под воздействием в виде сейсмической волны. Влияние волны классически воздействует на балку в местах ее крепления. В качестве функции, имитирующей поведение сейсмических волн, используется функция материнского вейвлета “mother wavelet”, имеющая вид:

${\Psi }_{a,b}\left(t\right)=\Psi \left(\frac{t-b}{a}\right),$ (62)

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

Начнем анализ полученных результатов со сравнения классической и “гистерезисной” балок с гистерезисным преобразователем Боука — Вена. Как видно из рис. 16 на начальных временных интервалах (с 0 по 3 отчет модельного времени) балки ведут себя одинаково, однако спустя 8 интервала временного отчета “гистерезисная” балка стаблизируется, в то время как классическая балка демонстрирует периодические колебания. Такая структура колебаний объясняется диссипирующими свойствами гистерезисного преобразователя (рис. 17), который трансформирует энергию колебаний сейсмической волны в иные формы.

Рис. 16. Слева — график материнской вейвлет-функции в зависимости от времени (вейвлет Добеши 12 порядка с параметрами a=1, b=4). Справа — отклонение дискретной точки балки с координатами L/2 от положения равновесия в зависимости от времени: синяя кривая – классическая модель балки, оранжевая кривая — гистерезисная балка с использованием преобразователя Прандтля — Ишлинского, красная кривая — гистерезисная балка с использованием преобразователя Боука — Вена

Fig. 16. Deflection graphs of the point on the beam with coordinates x=L/2 depending on the time under the influence of the seismic load presented as the Daubechies wavelets of the 12th order with parameters a=1, b=4: the orange curve — Prandtl-Ishlinskii model, the blue curve — the classical model

Рис. 17. Параметрические графики гистерезисных петель (время t как параметр): входно-выходная зависимость входа u_xx от выхода Γ[u_xx] в точке x=L/2: а — преобразователь Прандтля — Ишлинского, б — преобразователь Боука — Вена

Fig. 17. Parametric graphs of hysteresis loops (time t as a parameter): input-output dependence of the input u_xx on the output Γ[u_xx] at the point : a — the Prandtl — Ishlinskii operator, b — the Bouc — Wen operator 

Также из результатов моделирования, представленных на рис. 16, видно, что реация классической и “гистерезисной” балок различна. А именно в моментах воздействия сейсмической волны гистерезисная балка демонстрирует возмущения меньшой амплитуды (особенно четко это наблюдается для балки с использованием модели Прандтля — Ишлинского) в сравнении с классической балкой. После прекращения сейсмического воздействия обе балки стабилизируются в состояние равновесия, однако гистерезисная балка, формализуемая преобразователем Прандтля — Ишлинского, оказывается в положении равновесия, отличном от исходного (рис. 16, б), что согласуется с упруго-пластическими свойствами гистерезисного материала.

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

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

Заключение

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

Финансирование. Статья поддержана РНФ. Грант № 23-29-00696.

Информация о конфликте интересов: автор и рецензенты заявляют об отсутствии конфликта интересов.

^[1]  Здесь $u_{xxxx}=\frac{{\partial }^{4}u}{\partial x^4}$, $u_{tt}=\frac{{\partial }^{2}u}{\partial t^2}$.

^[2] Также известное как “strain hardening”, которое представляет собой упрочнение материала (например, металла или полимера) за счет пластической деформации, которое появляется из-за возникновения и движения дислокаций в кристаллической структуре материала

^[3] Кривизна равна нулю только в том случае, когда каждая точка балки неподвижна, однако как только балка подвергается деформациям, то есть некоторые бесконечно малые элементы балки начинают отклоняться от положения равновесия, кривизна становится ненулевой.

^[4] Очевидно, что производную необходимо взять только у слагаемого ${\Gamma }_{PI}[u_{xx}+\epsilon {\eta }_{xx}]$, поскольку  не зависит от $\epsilon$.

^[5] Вследствие использования леммы первое интегральное слагаемое будет равно нулю их-за разности $\eta (x,t)-\eta (x,t_0)$.

About the authors

Evgeny A. Karpov

Author for correspondence.
Email: believedream95@gmail.com
ORCID iD: 0000-0001-7032-4375

postgraduate student, student of the Department of Digital Technologies

References

Supplementary files

Supplementary Files

Action

1. JATS XML

Download

2. Fig. 1. Illustration of the main difference between Timoshenko and Euler — Bernoulli beam theories  

Download (72KB)

Indexing metadata

3. Fig. 2. Graphical representation of the Prandtl operator. E and P are the elastic and the plastic elements respectively

Download (16KB)

Indexing metadata

4. Fig. 3. Dependence of the input ε(t) on the output σ(t) for the Prandtl operator with the threshold value h and the domain Ω=Ω(Гр, E) presented as a band with the slope angle determined by the coefficient E

Download (34KB)

Indexing metadata

5. Fig. 4. Inner structure of the Prandtl — Ishlinskii operator with the threshold values h=h1, h2, ..., h∞, distributed using a distribution function and connected in parallel by the Prandtl operators Гр=Гр1, Гр2, ..., Гр∞

Download (59KB)

Indexing metadata

6. Fig. 5. Schema of parallel connection of the stop-operator models in the Prandtl — Ishlinskii model

Download (63KB)

Indexing metadata

7. Fig. 6. Bending of the cross section of the beam during the deformation

Download (97KB)

Indexing metadata

8. Fig. 7. Deformed part of the beam obtained by means of two parallel cross sections Ω(x) and Ω(x+Δx)

Download (172KB)

Indexing metadata

9. Fig. 8. Grid  with size M for N. Undefiined grid nodes are denoted by white color, nodes denoted by black color are calculating based on boundary conditions  

Download (91KB)

Indexing metadata

10. Fig. 9. Fixed boundary of the beam taking into account the external load g(x,t) (control) and seismological load S(t) at the boundary connection  

Download (50KB)

Indexing metadata

11. Fig. 10. Graphical representation of the load function

Download (291KB)

Indexing metadata

12. Fig. 11.Vibration graphs at the point x=L/2 depending on time: the blue curve denotes the classical beam, the orange curve denotes the Prandtl — Ishlinskii model, the red curve denotes the Bouc — Wen model ( а, b); parametric graphs of hysteresis loops (time t as a parameter): input-output dependence of the input uxx on the output Γ[uxx] at the point x=L/2 ( c, d) 

Download (671KB)

Indexing metadata

13. Fig. 12. Fourier spectrum of the vibrations of the beam at the point x=L/2: left — classical model, middle — Prandtl — Ishlinskii model, right — Bouc — Wen model 

Download (136KB)

Indexing metadata

14. Fig. 13. Vibration graphs at the point x=L/2 depending on time under the impulse load with the force A=1000 and parameters χ=0.5, τ=0.5: the blue curve denotes the classical beam, the orange curve denotes the Prandtl — Ishlinskii model, the red curve denotes the Bouc — Wen model 

Download (269KB)

Indexing metadata

15. Fig. 14. Parametric graphs of hysteresis loops (time t as a parameter): input-output dependence of the input uxx on the output Γ[uxx] at the point x=L/2: a — the Prandtl — Ishlinskii operator, b — the Bouc — Wen operator 

Download (204KB)

Indexing metadata

16. Fig. 15. The red curve corresponds to the position of the hysteretic beam formalised using the Bouc — Wen operator, the blue curve corresponds to the position of the classical beam, and the green curve corresponds to the point of impact  

Download (92KB)

Indexing metadata

17. Fig. 16. Deflection graphs of the point on the beam with coordinates x=L/2 depending on the time under the influence of the seismic load presented as the Daubechies wavelets of the 12th order with parameters a=1, b=4: the orange curve — Prandtl-Ishlinskii model, the blue curve — the classical model

Download (166KB)

Indexing metadata

18. Fig. 17. Parametric graphs of hysteresis loops (time t as a parameter): input-output dependence of the input uxx on the output Γ[uxx] at the point : a — the Prandtl — Ishlinskii operator, b — the Bouc — Wen operator 

Download (242KB)

Indexing metadata