CC BY
0
МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ
Б01: 10.18287/2541-7525-2024-30-3-35-62 УДК 510.6
Динамика балки Эйлера — Бернулли с учетом распределенных гистерезисных свойств
Карпов Е.А.
Воронежский государственный университет, г. Воронеж, Российская Федерация; [email protected];
Аннотация. В настоящей статье изучается математическая модель балки с распределенными гистерезисными свойствами. Свойства гистерезиса формализуются в рамках двух подходов: феноменологического (модель Боука — Вена) и конструктивного (модель Прандтля — Ишлинского). Уравнения колебаний балки получены с использованием известного подхода Гамильтона. Рассмотрены динамические характеристики балки с распределенным гистерезисом при различных видах внешней нагрузки: импульсной, периодической и сейсмической. Численное моделирование показывает, что балка гистерезиса более "устойчива" к внешним нагрузкам, чем классическая балка Эйлера — Бернулли. Эти результаты могут найти применение в области проектирования сейсмостойких конструкций и зданий.
Ключевые слова: гистерезис; балка Эйлера — Бернулли; балка; модель Боука — Вена; нелинейная динамика, устойчивость.
Поступила: 05.06.2024 Рассмотрена: 02.08.2024 Принята: 02.09.2024
Научная статья
Введение
В настоящее время проектирование несущих конструкций различных инженерных сооружений, а также разработка адекватных физических моделей являются неотъемлемой частью строительного процесса. Особенно остра задача выбора конструктивных составляющих зданий стоит в городах и странах, находящихся в сейсмоопасных регионах. Одним из важнейших этапов моделирования несущих конструкций в этом случае является идентификация форм их поведения под воздействием внешних возбуждений (нагрузок). Обычно классификация нагрузок выполняется по нескольким критериям: по характеру приложения (сосредоточенные и распределенные), по продолжительности во времени (переменные и постоянные), по характеру действия (статические и динамические). В зависимости от типа нагрузки проявляются те или иные свойства несущих конструкций, являющиеся, в свою очередь, следствием внутренних особенностей используемых материалов.
Рис. 1. Графическое представление принципиального различия между теориями Тимошенко
и Эйлера — Бернулли
Fig. 1. Illustration of the main difference between Timoshenko and Euler — Bernoulli beam theories
Известно, что одним из базовых элементов любой несущей конструкции является балка. Математическое описание балки обычно формализуется в рамках теории Эйлера — Бернулли [1] или теории Тимошенко [2] (рис. 1). Основное различие между вышеозначенными подходами заключается в поведении балки при поперечных деформациях. В теории Эйлера — Бернулли предполагается, что во время деформаций поперечное сечение перпендикулярно нейтральной оси. В то же время в теории Тимошенко оно не перпендикулярно нейтральной оси и напрямую зависит от интенсивности деформации. В настоящей статье модель колебаний балки строится с использованием теории Эйлера — Бернули в предположении, что балка имеет шарнирное закрепление в двух концах и обладает нулевыми начальными условиями (балка покоится и не имеет начальной скорости). Следуя работе [3], классическое уравнение колебаний балки имеет следующий вид:
щь + а2Чхххх = д(х, Ь), Ь > 0, х е [0, Ь],
и(х,0) = 0, щ (х,0) = 0, (1)
' и(0, Ь) = и(Ь, Ь) = Б(Ь),
Щхх (0, Ь) = ихх (Ь, Ь) = 0.
где и = и(х, Ь) - отклонение балки от положения равновесия, индексы обозначают частные производные1, д(х, Ь) есть функция, характеризующая вертикальную нагрузку (управление). Б(Ь) - функция, характеризующая внешнее воздействие, прикладываемое к закрепленным концам балки, а2 = Е1 - произведение модуля упругости и момента инерции, Ь > 0 есть время моделирования, Ь - длина балки.
Балки в зависимости от своих конструктивных особенностей находят свое применение в различных элементах строительных (и не только) конструкций. Например, двутавровые балки используются для обеспечения устойчивости перекрытий между этажами и в мостовых сооружениях, швеллерные же балки используются для армирования железобетонных конструкций, а также в качестве рам, каркасов и т. д. Кроме того, балки по-разному реагируют на внешнюю нагрузку, которая, в свою очередь, может быть импульсной, вибрационной, подвижной [4], ударной, ветровой. При этом основная задача, решаемая при выборе типа балки применительно к задачам строительства, заключается в обеспечении устойчивости всего сооружения по отношению к потенциальным внешним воздействиям, к числу которых, в первую очередь, относятся сейсмические волны как природного, так и техногенного характера. В таком случае балка выступает в роли демпфирующего элемента, трансформируя входящую энергию воздействия в иные ее формы [5].
Известно, что во время нагрузок балка начинает деформироваться, порождая во внутренней структуре материала силы, стремящиеся вернуть ей исходную форму. Если после снятия нагрузки она принимает форму, отличную от исходной, то следует говорить об упруго-пластическом гистерезисе. Он играет важную роль в современных исследованиях реальных систем и процессов. Это явление было обнаружено во многих областях современной науки: физике [6-14], экономике, биологии, химии и т. д. Недавние исследования показывают важность и целесообразность использования гистерезиса в нейронных сетях, например, включая его как звено или как самостоятельный элемент в функцию активации соответствующей сети. Такой широкий спектр применения обуславливается важными структурными особенностями гистерезиса - зависимостью текущего состояния системы от предыстории и начального состояния как от параметра. Еще более важной особенностью гистерезиса является тот факт, что в ряде задач он позволяет моделировать процесс диссипации энергии, что, в свою очередь, позволяет использовать гистерезисные преобразователи в качестве управляющего элемента, стабилизирующего поведение динамических систем [15-17].
Теоретические модели, описывающие гистерезисные явления, восходят к классическим работам Прандтля [18], Прейзаха [19] и др. В начале 60-х годов XX века появились первые феноменологические модели [20-23]. Теория конструктивных моделей гистерезиса, учитываю-
-I „ d4u d2u
Здесь uxxxx = dX4 , utt = •
щая физические особенности конструктивных элементов исследуемых систем, получила свое развитие в трудах М.А. Красносельского и его учеников в Воронежском университете в 60-х годах прошлого века. Фундаментальная монография [24] с подробным описанием свойств конструктивных моделей гистерезиса вышла в 1983 г. Разработанная в ней методология позволила описывать гистерезис на языке нелинейных операторов, зависящих от своего начального состояния как от параметра. Математический аппарат, предложенный в ней, базируется на выделении элементарных носителей гистерезисных свойств - гистеронов - преобразователей, наделенных пространством состояний и входно-выходными соответствиями. Такой подход позволил имплементировать модели гистерезиса в состав более сложных систем, описываемых дифференциальными уравнениями.
В настоящей статье предлагается модель балки с учетом распределенных гистерезисных свойств, в разд. 1 описываются используемые в работе модели гистерезиса, разд. 2 посвящен выводу уравнений колебаний гистерезисной балки, в разд. 3 описаны численные методы, используемые в работе, а в разд. 4 приводятся результаты моделирования и сравнительный анализ классической и гистерезисной балок. В качестве моделируемой нагрузки будут выступать периодические и импульсные воздействия, а также нагрузка сейсмического характера.
1. Модели гистерезиса: феноменологический и конструктивный подходы
1.1. Преобразователь Прандтля
Одной из самых известных и простых гистерезисных моделей является преобразователь Прандтля [25-27] (также известный как "stop operator"). Он используется для моделирования упруго-пластических волокон [28], где состояния полностью определяются величинами переменной деформации и переменного напряжения, а также применяется в качестве элементарной составляющей в сложных континуальных гистерезисных моделях в рамках теории упруго-пластического гистерезиса. Действительно, преобразователь Прандтля идеально описывает поведение материала: материал, находящийся в изначальном состоянии, когда нормальное напряжение и деформация отсутствуют, под нагрузкой начинает вести себя упруго до некоторого порогового значения h, а после демонстрирует пластические свойства, находясь под постоянным напряжением. Это поведение также называется "пластическое течение", так как при уменьшении напряжения материал обладает пластическими свойствами, пока не достигнет порогового значения с противоположным знаком —h. Однако отметим, что преобразователь Прандтля позволяет моделировать упруго-пластические свойства без деформационного упрочнения2.
Структурно этот преобразователь может быть представлен как последовательное соединение линейного упругого элемента E (например, пружина) с жестким, идеально пластическим элементом P (например, брусок, двигающийся по поверхности, под действием силы сухого трения) (рис. 2). Если предположить, что упругость пружины E равна единице, следует говорить о более простом преобразователе - упоре.
Описание преобразователя Прандтля начнем с определения его состояний. Состояниями этого преобразователя являются пары "вход-выход" {x, u} с областью определения в виде полосы П = П(Гр, E), угол наклона которой зависит от коэффициента E, определяющего модуль упругости. Причем полоса П зависит также и от порогового значения h (полоса пересекает ось Ox в точках ±h), которое является еще одной характеристикой этой модели. Иначе говоря, модель может быть формализована в виде равенства:
o(t) = Гр [h]e(t),
2Также известное как "strain hardening", которое представляет собой упрочнение материала (например, металла или полимера) за счет пластической деформации, которое появляется из-за возникновения и движения дислокаций в кристаллической структуре материала
-ЛМАА—s—
Рис. 2. Графическое представление гистерезисного преобразователя Прандтля. Здесь E и P — упругий
и пластический элементы соответственно
Fig. 2. Graphical representation of the Prandtl operator. E and P are the elastic and the plastic elements
respectively
где
ГР[%(£) = <(*' Т)Ь если е(Ь) > е(Ь - т),
1тах{-й, р(Ь, т)}, иначе
р(Ь, т) = Е [е(Ь) - е(Ь - т)] + о(Ь - т) или в более компактной форме:
о(Ь) = Гр [Н]е(Ь) = тт{й,тах{-й, р(Ь, т)}}, (2)
где о(Ь - т) — значение выхода гистерезисного преобразователя в момент времени Ь - т, е(Ь) -входной сигнал, о(Ь) — выходной сигнал (реакция преобразователя). Графическое представление описанного выше преобразователя отображено на рис. 3.
Рис. 3. График зависимости входа e(t) от выхода o(t) для преобразователя Прандтля с пороговым значением, равным h, и областью определения П = П(Гр, E) в виде полосы с углом наклона,
зависящим от коэффициента E
Fig. 3. Dependence of the input e(t) on the output o(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). В настоящем разделе приводится математическое описание этого преобразователя.
г; г:
ЛДМЛ—■
■ЛМААЛМ-
■ЛМДА-
h, h2 h,
Рис. 4. Внутренняя структура континуального преобразователя Прандтля — Ишлинского с пороговыми значениями h = hi, h2, •••,hTC, распределенными с помощью некоторой функции распределения и параллельно соединенными преобразователями Прандтля Гр = Гр1, Гр2, •••, Гртс
Fig. 4. Inner structure of the Prandtl — Ishlinskii operator with the threshold values h = hi,h2, •••,hTC, distributed using a distribution function and connected in parallel by the Prandtl operators
ГР = ГР^ •••, rP~
Зададим семейство однопараметрических преобразователей Прандтля (2), зависящих от порогового значения h как от параметра. Пусть параметр h, характеризующий пороговое значение, принадлежит некоторому множеству Ф. Пусть задана неубывающая непрерывная слева функция S = S(h), h > 0 на множестве Ф, удовлетворяющая следующим двум условиям:
lim E(h) = 0,
Г
Jo
\E(h)\dh < œ.
Также пусть задано множество Z непрерывных функций о (И), удовлетворяющих условию
\о(И) | <И.
Теперь определим множество П(Гр, Е) возможных состояний преобразователя Гр как
П(Гр, Е) = {{£, 01 (И)}, (е, 02 (И)},..., (е, о^(И)}},
где е — это вход, ог(И) — выход (реакция) г-го преобразователя Грг (каждый преобразователь зависит от порогового значения как от параметра). Входно-выходные соответствия такого преобразователя описываются следующим равенством:
o(t) = rpi [to, oo(h); Z]£(t) = Гр [to, ffo(h); h]£(t)dS(h),
Jo
а переменное состояние {£(t),a(h,t)} определяется равенством
о (h; t) = Q[to, zo(h); S]£(t) = Гр [to, oo(h); h]£(t),
(3)
(4)
где Q есть пучок гистеронов (конечный или в общем случае бесконечный набор гистеронов), преобразователь Гр описывает параллельное соединение бесконечного набора преобразователей,
описываемых соотношением (2) (см. рис. 4), оо(к) стоит понимать как значение выхода в предыдущий момент времени (альтернативная запись — т,к)).
Согласно определениям и свойствам, изложенным в работе [24], подобный преобразователь обладает свойствами статичности и виброкорректности, однако свойствами управляемости и детерминированности не располагает. В ряде задач, вместо континуального преобразователя
z,(h)
Рис. 5. Блок-схема параллельного соединения упоров в модели Прандтля — Ишлинского
Fig. 5. Schema of parallel connection of the stop-operator models in the Prandtl — Ishlinskii model
Прандтля — Ишлинского удобно использовать его конечномерный аналог (именно такая ситуация имеет место в компьютерном моделировании гистерезисных зависимостей). Более того, в моделях реальных физических систем практически всегда реализуется конечномерный аналог континуального преобразователя. Пусть задан конечный набор размера N £ Z гистеронов Гр, формализуемых посредством преобразователя Прандтля. Предполагая, что гистероны соединены параллельно (как и в континуальной модели), и переходя от интеграла к интегральной сумме, преобразователь примет вид (рис. 5):
N
°(t) = Yj f (hi )Гр [to, a; hi]e(t). (5)
i=1
Большое разнообразие работ, использующих преобразователь Прандтля — Ишлинского, можно условно разделить на две категории: фундаментальные и прикладные. В рамках фундаментальных работ исследуются структурные особенности преобразователя [29-31], вопросы корректности, а также проблемы идентификации его параметров [32]. В работах прикладной направленности преобразователь используется для описания гистерезисных явлений в механических системах, например актуаторы, магнитореологические эластомеры [33], несущие конструкции [34] и т. д.
1.3. Феноменологический подход: преобразователь Боука — Вена
Наряду с конструктивным, не менее известным является феноменологический подход к моделированию гистерезиса. Он позволяет моделировать петли гистерезиса разной структуры, не принимая во внимание физическую основу моделируемого процесса. Обычно феноменологический подход применяется в том случае, когда необходимо оценить реакцию системы на сигнал лишь на основе экспериментальных данных о характеристиках входного воздействия и отклике системы на него. Одним из самых известных преобразователей в рамках этого подхода является преобразователь Боука — Вена, предложенный в [21; 22]. Входные-выходные соответствия преобразователя Боука — Вена формализуются уравнением:
Гbw = x A - [в sign (XfBw) + y] |rBwГ},
(6)
где х = х(Ь) — входной сигнал, Гвм = Гвм— выходной сигнал (реакция преобразователя), у, п, в, А > 0 — безразмерные параметры. Как можно заметить, преобразователь Боука — Вена является четырехпараметрическим, что позволяет моделировать различные по форме и площади гистерезисные петли. Отметим, что проблеме идентификации этих параметров посвящено достаточно большое количество публикаций.
Очевидным преимуществом преобразователя Боука — Вена являются простота численной реализации и развитый математический аппарат [35-38]. В этих работах рассмотрены задачи возникновения предельного цикла, взаимная зависимость параметров преобразователя, адаптация модели применительно к некоторому набору прикладных задач и т. д. На сегодняшний день преобразователь Боука — Вена применяется и для моделирования процессов различного масштаба [15; 39].
2. Уравнение колебаний гистерезисной балки: гамильтоновский подход
В настоящем разделе производится вывод уравнений колебаний балки с распределенными гистерезисными свойствами (рис. 6). Отметим, что в подавляющем большинстве прикладных задач рассматривались объекты с сосредоточенными гистерезисными свойствами. Распределенные гистерезисные свойства рассматривались в весьма ограниченном количестве работ, из которых отметим работу [34]. Ниже приводится вывод уравнений колебаний гистерезисной балки, основанный на вариационном принципе, а именно принципе наименьшего действия Гамильтона. В рамках этого подхода искомые уравнения являются следствием экстремального значения функционала действия.
Рассмотрим стержень, работающий на изгиб, имеющий постоянную площадь поперечного сечения S = const, симметричную относительно обеих осей Oy и Oz, при этом нейтральная ось стержня совпадает с осью Ox. Отметим, что в текущем выводе учитываются только поперечные колебания в плоскости Oxy, продольными колебаниями пренебрегают. Пусть отклонение гисте-
Положение Положение
поперечного сечения покоящейся балки
деформации
Рис. 6. Изгиб поперечного сечения балки во время деформации Fig. 6. Bending of the cross section of the beam during the deformation
резисной балки от положения равновесия описывается следующим функционалом действия:
I(u) = f1 Ldt, (7)
Jt о
где to < t < ti — некоторый промежуток времени, u = u(x, t) — кривая, описывающая динамику
балки, L — функция Лагранжа, равная разности кинетической и потенциальной энергии:
L = T — и, (8)
где U — потенциальная энергия, T — кинетическая энергия. Определим потенциальную энергию балки как ^
U = II Mxdx, (9)
2 «70
где M — изгибающий момент, а х — деформация балки в точке x. Деформацию определим через кривизну.
Пусть кривая, описывающая положение балки в плоскости OXy, задана уравнением F = = F(x, у) = 0 (— £ К2), тогда ее кривизна будет описываться уравнением:
\—2—хх — FxFyFxy + —Х Fyyl к — _
— + -2у )3/2 •
Так как продольными колебаниями можно пренебречь (согласно предположению), то уравнение, описывающее положение балки, будет иметь вид —(х,у) = у — и(х) = 0, а ее кривизна определяться следующим соотношением:
(1 - «J)3/2'
Ограничиваясь рассмотрением лишь малых углов отклонения поперечного сечения ф ^ 1 (ux ^ 1 и как следствие1 - uX —* 1), кривизна будет определяться равенством:
К = uxx-
В свою очередь, изменение кривизны балки равно3
X = k - k0 = k - 0 = k = uxx (10)
где ko — начальная кривизна покоящейся балки, равная нулю.
Для вывода уравнения изгибающего момента выберем два параллельных сечения П(х) и П(х + Дх) на расстоянии Дх. Очевидно, что углы поворота в этих сечениях будут различными во время деформации. Вследствие этого происходят растяжение (в верхней части балки) и сжатие (в нижней части балки) материала балки в продольном направлении. Рассмотрим часть балки, находящейся на расстоянии r = r(y) от ее оси симметрии. Как видно на рис. 7, длина растянутого материала изменилась на Дв = гДф. Относительное удлинение материала
Дв Дф
е(х, Дх, y, t) = - — = -r—1, (11)
к J ' Дх Дх
а деформация растяжения-сжатия в точке х
г(х, y, t) = - lim г(х, Дх, y, t) = -r lim —Ф = -гфх. (12)
v V ' Дх^0 у Дх^0 Дх т
Угол отклонения поперечного сечения ф от "нормального" положения в точке х (см. рис. 6) с учетом ортогональности оси симметрии направлению деформации равен:
Ф = UX. (13)
3Кривизна равна нулю только в том случае, когда каждая точка балки неподвижна, однако как только балка подвергается деформациям, то есть некоторые бесконечно малые элементы балки начинают отклоняться от положения равновесия, кривизна становится ненулевой.
«
XX
К
Рис. 7. Деформированная часть балки, полученная двумя параллельными сечениями П(х)
и П(х + Ax)
Fig. 7. Deformed part of the beam obtained by means of two parallel cross sections Q(x) and Q(x + Ax)
Таким образом, подставляя (13) в (14), получим
е(х, у, Ь) = -ruxx■ (14)
Изгибающий момент, возникающий относительно оси Oz вследствие нормального напряжения в сечении П(х), описывается уравнением
M(x, t) = radS. (15)
Ja
Теперь вычислим нормальное напряжение в сечении П(х), полагая, что материал балки обладает упруго-пластическими свойствами с деформационным упрочнением. Для этого воспользуемся преобразователем Прандтля — Ишлинского:
✓-»то
ст(Ь) = Ги^о,сто;Е]е(ь) := Гр[Ьо,сто;И.]е(^Е(К), (16)
где ст(Ь) — выход (реакция преобразователя), ГpI — преобразователь Прандтля — Ишлинского, е(Ь) — вход. Выполним переход от одномерного преобразователя Прандтля — Ишлинского к двумерному, введя пространственный параметр С, принадлежащий множеству измеримых параметров, по Лебегу, как предлагается в работе [34]. Тогда параметрический преобразователь примет вид:
✓-»то
ст(С,Ь) = Гр1 [Ьо,сто;Е]е(С,Ь) := Гр[Ьо,сто;Н]е(С,К),
о
где С — параметр, характеризующий пространственную (по материалу балки) плотность распределения элементарных носителей гистерезиса. Или, используя следующее соотношение С = (х,у) (параметр С представлен в виде пространственной переменной), окончательно получим:
✓-»то
ст(х,у,Ь) = Гр1 [Ьо,сто;Е]е(х,у,Ь) := ГР[Ьо,сто;Н]е(х,у,Ь^Е(х,у,К). (17)
о
С учетом соотношения (17) уравнение (15) примет вид:
М(х,Ь) = ^ г^ Гр[Ьо,сто;К]е(х,у,Ь^Е(х,у,К)^dS. (18)
Раскрыв скобки и учитывая соотношение (14), получим
И.
гГр [Ь0, оо; Н] (-тихх)йЕ(х, у, Н)йБ. (19)
Теперь воспользуемся свойством преобразователя Прандтля [34], а именно соотношением
тГр[Н]{-тихх) = -т2Г[Н/\т\]ихх. (20)
Для подынтегрального выражения из соотношения (19) с учетом соотношения (14)
М(х, Ь) = ^ -,(Х Гр[£о, оо;Н/\т\]иххйЕ(х,у,Н)^йБ (21)
и введя переменную Е(х,Н) = ^ т2Е(х, у, \ т\Н)йБ, получим
М(х, Ь)= - Гр[£о, оо; Н]иххйнЕ(х, Н). (22)
о
Как можно заметить, уравнение (22) совпадает с уравнением (16), что позволяет записать соотношение в следующей форме:
М(х, Ь) = -Гиихх. (23)
Таким образом, связывая отношение деформации к напряжению посредством преобразователя Прандтля — Ишлинского (17), по сути, осуществляется распределение гистерезисных свойств вдоль длины балки.
Окончательно принимая во внимание соотношения (10) и (23), уравнение для потенциальной энергии примет вид
и = -1 ( иххГр\[ихх]йх. (24)
2 «70
Кинетическую энергию определим классическим образом:
If p(x)Si
T = p(x)Sufdx. (25)
Тогда интеграл действия для такой системы примет вид
"*2 1 / г1 г1
(,==1 ^f uxxГР1 [uxx]dx + f pSufdx^dt.
(и) = I 2\ ] иххги[ихх]йх + I рБи2йх]йЬ. (26)
Из принципа наименьшего действия в форме Гамильтона следует, что функция и = и(х, Ь) является экстремалью функционала (26). Рассмотрим вариацию этой функции и(х, Ь) = и(х, Ь)+ +ец(х,Ь), для которой справедливо следующее:
• в начальный и конечный моменты времени функция и(х, Ь) не варьируется;
• ц(х, Ь\) = ц(х, Ь2) = 0 для таких х, для которых верно 0 < х < Ь. Таким образом, функционал (26) для кривой Еи примет вид
I = ^ ^^ (ихх + ецхх)Г[ихх + ецхх]йх + рБ ^ (щ + ещ )2йх^йЬ. (27)
Ясно, что функционал действия (27) достигает экстремума в точке е = 0. Согласно необходимому условию существования экстремума, положим вариацию функционала действия относительно е равной нулю
йт
0. (28)
dl_ de
e=0
Перед вычислением производной по £ раскроем скобки внутри подынтегральных уравнений (27):
Uxx Грт [uxx + ецхх ] + еПххГт [uxx + Щхх ]dx,
J 0
(29)
0
pS f (u2 + leutrjt + n2)dx. (30)
0
Подставляя в (28) выражение интеграла действия (27) с учетом (29) и (30), вычислим производную по £:
Мг
АЬ
чч- —/ I -/С
где
Г2 ((Ii + h) + pSh)edt = 0, (31)
Л
Ii = U xx Гр1 [u xx + eqxx]dx, (32)
0
о
1г = £ЦххГр1 [и хх + £Цхх ^х, (33)
о
1з = Г (и2 + 2£ЩЩ + ¥х. (34)
о
Предположим, что функция ц(х, Ь) на концах равняется нулю, то есть п(о, Ь) = ц(Ь, Ь) = о. Рассмотрим более детально каждое слагаемое (32)-(34).
Найдем частную производную по £ для (32). Для этого продифференцируем подынтеграль-
4
ное слагаемое :
uxx (Грт [u xx + enxx])e dx. (35)
0
Для дальнейших вычислений предположим, что весовая функция Е дифференцируема и ее производная есть функция f (h). Следуя этому, преобразователь Прандтля — Ишлинского примет вид:
✓-»то
ГРТ [h]e = f (h)rP [h]edh. (36)
0
Учитывая (36), запишем полную форму производной преобразователя Прандтля — Ишлинского из соотношения (35):
✓-»то
(Грт [uxx + en xx ])e = f (h)(rp[uxx + enxx ])e dh, (37)
0
где Гр — преобразователь Прандтля. Представляя преобразователь Прандтля в форме (2):
ГР [<p(x, t, e)] = min{h,max{-h, x, t, t0, e)}}, (38)
где x, t, t0,e) = Гр0 + §(x, t,e) - §(x, t0,e) и §(x, t,e) = uxx + er\xx, его производная примет вид:
ru,, , i0, если t, to, e) e [-h Щ,
Гр ^ e)]e = ^ . - . - (39)
[nxx(x, t) - nxx(x, t0), иначе. Полученный выше результат можно представить в виде
Гр[<p(x, t,e)]e = (nxx(x, t) - nxx(x, t0))x
x (&(^(x, t, t0,e) + h) - ©($(x, t, t0,e) - h)), (40)
4Очевидно, что производную необходимо взять только у слагаемого Гр1 [ихх + £Цхх], поскольку и не зависит от £.
где ©(•) — функция Хевисайда. Тогда, подставляя полученный результат в (37) и вынося слагаемые, не зависящие от Н, получим
(ГИ [и хх + ец хх ])е =(П хх (х, Ь) Цхх (х/ Ь0 ))х
✓"»ТО
х /(Н)(@(у(х, Ь, Ь0, е) + Н) - ®(у(х, Ь, Ь0, е) - Н))йН. (41) Подставив в (35) выражение (41), обозначив интеграл как некоторую функцию в(х, Ь, Ь0, е):
ихх (Цхх
(х, Ь)
Цхх (х/ Ь0))в(х, Ь, Ь0, е)йх, (42)
ю
проинтегрируем два раза по частям, принимая за V функцию от ц (чтобы избавиться от двух производных по х), и окончательно получим
ЛЬ
(х,Ь) - п(х,Ь0))а(х,
f
0
(h)e = (n(x, t) - n(x, t0))a(x, t, t0,e)dx, (43)
0
где а(х, Ь, Ь0, е) — некоторая произвольная функция, полученная в ходе интегрирования. Теперь рассмотрим второй интеграл 12. Дифференцируя интеграл по е:
nxx Гр1 [uxx + enxx ]dx + e^xx (Гр1 [uxx + e^xx])edx (44)
00
f
0
и полагая £ = 0, получим следующее выражение:
nxxrpi[ixx]dx. (45)
0
Теперь продифференцируем выражение два раза по частям и окончательно получим
(I2 )e = nrpi [ixx ]xxdx. (46)
0
Нетрудно заметить, что в интеграле (34) только одно слагаемое зависит от е. Следовательно, дифференцируя по е, имеем:
2иЬг[Ьйх. (47)
Ю
Теперь, продифференцировав один раз по частям полученное выражение, окончательно имеем:
Тз)е =
Г
0
(Is)e = f 2quttdx. (48)
0
Таким образом, подставляя в (31) выражения (43), (46), (48), получим: J (J" (n(x,t) - n(x,t0))a(x,t,t0,0)dx + vp[uxx]xxdxj +
+ pSuttr[dxdt = 0. (49)
Заметим, что вышеполученное уравнение есть необходимое условие Эйлера локального экстремума функционала действия. Воспользуемся леммой Дебуа-Реймона, которая является основой для построения теории обобщенных функций, применительно к уравнению (49) и получим окончательно уравнение колебаний гистерезисной балки5:
Г[ихх]хх + рБиЬЬ = 0. (50)
5Вследствие использования леммы первое интегральное слагаемое будет равно нулю их-за разности П(х, Ь) - п(х, Ь0).
Учитывая потенциал внешней нагрузки и полагая произведение распределением массы на площадь поперечного сечения, равного единице, перепишем уравнение в виде
Щь + Г[ихх]хх = д(х, Ь), (51)
где д(х, Ь) — функция внешней нагрузки.
3. Численные методы: явная разностная схема и метод Рунге — Кутты
Задача нахождения решения дифференциального уравнения (в общем случае системы дифференциальных уравнений) является самой первой задачей в рамках исследования любой системы, описываемой дифференциальными уравнениями. Как правило, обычно используются точные аналитические методы решения. Ярким примером служит применение метода разделения переменных Фурье в задаче о колебаниях струны. Преимущество точных методов состоит в том, что они позволяют получить решение уравнения в виде некоторой комбинации элементарных функций (иногда решение представляется в виде квадратур от элементарных функций). Однако большинство практических задач решить с использованием этих методов не является возможным. Вызвано это либо особенностями области, в которой ищется решение уравнения (многие практические задачи накладывают ограничения на пространство параметров и области решений), либо вообще невозможностью нахождения аналитического решения. В таком случае удобно использовать численные методы, позволяющие искать решение, аппроксимируя исходное уравнение. Математический аппарат численных методов достаточно развит для его применения во многих задачах. Однако применение этих методов становится нетривиальной задачей, когда уравнения содержит операторные нелинейности гистерезисного типа. Для решения таких уравнений необходимо модифицировать классические методы.
В рамках настоящей статьи рассматриваются численные методы, преобразующие дифференциальное уравнение в частных производных с операторной гистерезисной нелинейностью к конечной системе алгебраических уравнений. Наиболее часто в литературе используется следующая интерпретация разностных схем: явные и неявные разностные схемы. Неявная разностная схема предполагает преобразование исходного уравнения к системе алгебраических уравнений с трехдиагональной матрицей с последующим решением (например, методом прогонки). Явная схема более проста в реализации и позволяет вычислять значение функции. опираясь на данные, полученные на предыдущем шаге. Неявные схемы, как правило, более устойчивы, чем явные. Однако, полагая отношение шага по времени много больше шага по пространству, можно добиться необходимой устойчивости. Ввиду сложности формализации гистерезисных операторов и их множественной суперпозиции в неявной разностной схеме, в настоящей работе используется явная разностная схема.
Введем в области П = {0 < х < L, 0 < Ь < Т} равномерную сетку с шагом по пространству Кх и с шагом по времени Кь, равную
Т T
Кх = N К = М (52)
где М — количество узлов сетки по временной координате, Т — время моделирования, N -количество узлов сетки по пространственной координате (рис. 8). Причем шаги подчиняются соотношению К ^ Кх для устойчивости описываемой схемы. Первые два временных шага Ь = о, Ь = 1 вычисляются на основе начальных условий. Используя граничные условия, нетрудно идентифицировать все значения в узлах с координатами (х(, Ь) где I = о, 1, N - 1, N - 2, ] = о, М. Для вычисления остальных значений аппроксимируем уравнение колебаний гистерезисной балки.
Уравнения, аппроксимирующие вторые производные по времени и пространству, примут вид:
Щ-1 - гЩ + Щ+1
ии = ->--2-—, (53)
M M-l M-2
9 •
О О О О О О
о о о о о о
о о о о о о
о о о о о о
о о о о о о
о о о о о о
о о о о о о
о о о о о о
о о о о о о
• •—• •
0 12 •
•-•-•-• X
N-2 N-1 N
Рис. 8. Область П = {0 < x < L, 0 < t < T} с введенной равномерной сеткой размера M на N. Узлы сетки, обозначенные белым цветом, вычисляются в ходе решения уравнения, узлы, обозначенные черным цветом, вычисляются из граничных условий
Fig. 8. Grid П = {0 < x < L,0 < t < 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
uj , - 2uj + u1.,, i-i i i+i ,Г.Л
uxx =-h2-' (54)
nx
где u^ соответствует значению в области П с координатами (i, j). Аппроксимация нулевых начальных и граничных условий очевидна, поэтому не приводится в настоящей статье. Аппроксимация функции нагрузки g(x, t) выполнена следующим образом:
g(x, t) = g( jhx, iht), (55)
а аппроксимация гистерезисного оператора имеет вид:
r[uxx] 1 - 2Г[uxx] + r[uxx]i+1
г [uxx]xx = -1--1-, (56)
hx
_j ul_ ,-2ul+ul+.
и, вводя новую переменную и. = —-г2——, согласно (54), окночательно получим соотношение:
1 hx
ги]i-i - 2Г[UUj]j + r[U]1+-
Г[„ i L 'H L 'H L 'H /гь7\ r[uxx] =--2-, (57)
hx
где Г — гистерезисный оператор. Подставляя полученные аппроксимации в уравнение (51), получим следующее соотношение:
ui-i - 2u;' + ui+i Г U]i-i - 2Г [U1 ]j + r[UU] i+i
u--— + —-hr-— = 9( jhx, 'ht). (58)
ht2 h2x
Оставляя в левой части уравнения только неизвестные, окончательно получим уравнение колебаний гистерезисной балки, аппроксимированное явной численной схемой:
uj+i = 2uj + -u^gdhx, 'ht) - h2 T[U' - + T[U' . (59)
hx
Таким образом, значение на j + i временном слое будет высчитываться на основе j и j - i слоев.
В зависимости от формализации гистерезисного звена будут рассматриваться различные подходы к его разностностной дискретизации. Так, для расчета реакции преобразователя Прандтля — Ишлинского будет использоваться его дискретный аналог (5). Выход преобразователя Боука — Вена рассчитывается, используя модификацию хорошо известного метода Рунге — Кутты четвертого порядка. Кратко опишем используемый алгоритм.
Пусть ¥ = ¥(£, ц) есть правая часть дифференциального уравнения преобразователя Боука — Вена со входом £ и значением этой функции в предыдущий момент времени ц. Введем вектор
п4, n2)
F(¿N-2, nN-2)
(60)
j _ ui-i-2ui +uj+i
j w+i и n ■ = Г 11 соответственно.
преобразователей Боука — Вена со значениями — 2 " >ц ~ ¿+1 Тогда модель Рунге — Кутты 4-го порядка применительно к решению уравнения Боука — Вена можно записать в векторном виде:
где
1
Г j+1 = rj + - (kij + 2k2j + 2кзj + k4j), 6
kij = htFj\i=i, Si Si
k j = hj={j+kij/2, j={j+k2j/ SH! +k3j '
(61)
k3j = htF\j={'+k2'/2, ^l 2l
k4j = htF j\
где Fj, kij, k2j, k3j, k4j есть векторы-столбцы на j временном интервале. Таким образом, представленный метод позволяет оптимизировать вычисления, значительно сократив время их выполнения (в средах, где есть поддержка векторного исчисления, например, среда Mathlab или библиотека NumPy для языка Python и др.).
4. Результаты моделирования колебаний
В настоящем разделе описываются результаты моделирования колебаний классической и гистерезисной балок в зависимости от разных типов нагрузок. Моделирование феномена гистерезиса осуществляется посредством двух преобразователей: Боука — Вена и Прандтля — Ишлинского. Также в разделе приводится сравнительная характеристика между этими моделями. Начальные условия полагаются нулевыми (балка покоится в начальный момент времени и не имеет начальной скорости соответственно) и ненулевыми граничные условия, когда моделируется сейсмическая нагрузка (на концах моделируется воздействие сейсмической волны) и нулевыми для иных нагрузок (рис. 9). Характеристики балки постоянны во всех экспериментах (кроме тех случаев, где явно указаны другие параметры):
• длина балки балки равна L = 1;
• время моделирования равно Т = 5;
• коэффициент а2 равняется о.5;
параметры численного моделирования полагались следующими:
• шаг по времени равен К = о.оо125;
• шаг по пространсву равен hx = 0.05; а параметры гистерезисных моделей:
• пороговые значения hi преобразователя Прандтля — Ишлинского распределялись в интервале (0,100) с шагом 0.5;
• коэффициент E преобразователя Прандтля — Ишлинского полагается равным 0.9.
• параметры преобразователя Боука — Вена полагаются равными: A = 0.9, ß = 0.2, у = 0.5, n = 0.2.
Отметим, что результаты были получены с использованием численных методов, описанных в разд. 3. Начнем описание полученных результатов с моделирования воздействия периодической структуры.
Load (control): g(x, t) ^ ^ ^ ^ ^ ^ ^ ^ ^ ^
A A
T T
S(t) Seismic influence: S(t) S(t)
Рис. 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
4.1. Моделирование нагрузки как функции периодической структуры
Одним из классических типов нагрузки является периодическая нагрузка. С ее помощью можно понять, как реагирует система на периодические воздействия с меняющейся частотой или амплитудой. Обычно в качестве периодической нагрузки принято использовать тригонометрические функции, изменяющиеся в зависимости от времени, вида A sin(ivt + фо) или A cos(Ш + фо), где A — амплитуда, ш — частота и фо — начальная фаза колебаний. В настоящей статье в качестве функции периодического воздействия будет использоваться функция g(x,t) = 0.12 sin(2n/2 - 4nt) (рис. 10).
Как можно видеть из полученных результатов (рис. 11), реакция каждого типа балки на периодическое воздействие разная. Классическая балка с начального момента времени начинает демонстрировать квазипериодические формы колебаний c двумя частотами (основной с наибольшей мощностью, и второстепенная, с меньшей), что подтверждается анализом спектра Фурье-преобразования (рис. 12, а). Балка с гистерезисными свойствами, формализуемыми преобразователем Боука — Вена, демонстрирует схожие результаты с классической балкой. А именно на начальном временном отрезке можно видеть схожие колебания с классической моделью. Однако, начиная с некоторого момента времени, устанавливаются периодические колебания с постоянной частотой (отметим, что частота у каждого из преобразователей отличается, ровно как и амплитуда колебаний). Анализ спектра Фурье-преобразования показывает, что использование преобразователя Боука — Вена способствует уменьшению мощности неосновной частоты колебаний (рис. 12, в) практически нивелируя ее. Используя в качестве носителя гистерезисных свойств преобразователь Прандтля — Ишлинского, можно добиться лучших
Рис. 10. График функции нагрузки g(x, t) = 0.12sin(n/2 - 4nt) Fig. 10. Graphical representation of the load function g(x, t) = 0.12sin(n/2 - 4nt)
M 0.000
J
- 0.005
a
- 0.010
-0.003-0.002-0.0010.000 0.001 0.002 0.003 rUxx(L/2, t)
в
-0.02 0.00 0.02 Г|„,и](Ь/2, t)
г
Рис. 11. Графики колебаний балки в точке x = L/2 в зависимости от времени: синяя кривая — классическая модель балки, оранжевая кривая — гистерезисная балка с использованием преобразователя Прандтля — Ишлинского, красная кривая — гистерезисная балка с использованием преобразователя Боука — Вена (а, б); параметрические графики гистерезисных петель (время t как параметр): входно-выходная зависимость входа uxx от выхода r[uxx] в точке 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 uxx on the output r[uxx] at the point x = L/2 (c, d)
а
результатов в сравнении с остальными моделями. А именно с ее помощью, во-первых, можно добиться уменьшения амплитуды основной частоты, во-вторых, редуцировать второстепенную частоту (рис. 12, б). Представленные гистерезисные зависимости (рис. 11, в, г) позволяют
оценить энергетический вклад гистерезисного звена в редуцирование колебаний: площадь петли пропорциональна энергии, "отводимой" гистерезисным звеном из системы.
Таким образом, преобразователь Прандтля — Ишлинского позволяет избавить колебания балки от квазипериодической структуры, трансформируя их в чисто периодическую, а преобразователь Боука — Вена способствует уменьшению мощности неосновной частоты колебаний балки.
0.025 .............................:
0.020
0.000 = ! ■.........................J
О 50 100 150 200 250 300 Fourier Spectrum
а
0.0005 .............................:
0.0004 0.0003 0.0002 0.0001
0.0000 .............................-
О 50 100 150 200 250 300
Fourier Spectrum
0.008 0.006 0.004 0.002
0.000 V.............................-
О 50 100 150 200 250 300
Fourier Spectrum в
Рис. 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. Моделирование воздействия импульсной нагрузки
Нередко в рамках прикладных задач важно знать реакцию несущей конструкции на импульсную нагрузку. Самым простым примером такой нагрузки является точечный удар по конструкции с некоторой "силой". Классическим методом формализации импульсной нагрузки является дельта-функция Дирака. Однако, принимая во внимание использование численных методов в работе, имеет смысл рассмотреть подобную ей функцию, которая в одномерном случае имеет вид:
х |Л, если Ь = 0,
<5(0 = I
10, иначе,
а в двумерном:
х I Л, если Ь = 0 и х = 0, 5(х, Ь) = \ '
10, иначе,
где Л — "сила" удара. Использование такой функции позволяет нивелировать бесконечности в расчетах, что сильно упрощает задачу моделирования. Ясно, что балка за счет внутренних сил должна "поглощать" удар, преобразуя энергию удара в иные ее формы, например, в тепловую энергею (нагрев балки). В настоящей статье исследуется реакция классической и "гистерезисной" балок на импульсную нагрузку (удар), приложенную к ее определенной дискретной точке в некоторый момент времени. Таким образом, функцию д(х, Ь) можно представить в виде:
д(х, Ь) = 5(х - х, Ь - т) = 5(х - х)<(Ь - т),
где х £ [0, Е] есть дискретная точка балки, в которой действует импульсная нагрузка, а т £ (0, Т) — момент времени, в который действует нагрузка.
Результаты моделирования, продемонстрированные на рис. 13, показали, что учет гистерезисных свойств приводит к совершенно иной динамике колебаний, нежели в классической модели. А именно начиная с момента времени Ь = т классическая балка начинает совершать колебания с постоянной частотой, причем в ходе экспериментов было установлено, что импульсная нагрузка любой "силы" выводит балку из состояния равновесия, заставляя балку
t t t
а б в
Рис. 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
бесконечно долго колебаться. Напротив, балка с учетом гистерезисных свойств с момента времени £ = т начинает отклоняться от положения равновесия, достигает некоторой точки, совершая при этом одиночное колебание малой амплитуды, после чего достигает устойчивого положения, отличного от положения равновесия и остается в этом положении бесконечно долго. При этом на рис. 13, в можно видеть, что результаты, демонстрируемые с использованием преобразователей Прандтля — Ишлинского и Боука — Вена, отличаются, однако схожи по своей структуре. Схожесть результатов подтверждается и подобностью гистерезисных петель, продемонстрированных на рис. 14.
- 0.10 - 0.05 0.00 - 0.12 - 0.10 - 0.08 - 0.06 - 0.04 - 0.02 0.00
rUxx(L/2, t) rUxx(L/2, t)
аб
Рис. 14. Параметрические графики гистерезисных петель (время t как параметр): входно-выходная зависимость входа uxx от выхода r[uxx] в точке x = L/2: а — преобразователь Прандтля — Ишлинского,
б — преобразователь Боука — Вена Fig. 14. Parametric graphs of hysteresis loops (time t as a parameter): input-output dependence of the input uxx on the output r[uxx] at the point x = L/2: a — the Prandtl — Ishlinskii operator, b — the Bouc — Wen
operator
Наблюдаемое поведение колебаний у балки с учетом гистерезисных свойств можно интерпретировать как более удачное, чем у балки, моделируемой классическим методом. Такой вывод можно сделать, полагая, что, получая удар, балка деформируется в сторону удара (соседние точки деформируются, в том числе образуя адекватный прогиб, см. рис. 15), при этом гистерезис способствует диссипации энергии удара, нивелируя дальнейшие колебания балки, в отличие от классической модели, колебания которой имеют периодический характер.
0.010
0.008
?ï 0.006 H
§ 0.004 0.002 0.000
Рис. 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", имеющая вид:
4a,b(t)= , (62)
где а и b являются параметрами растяжения и сдвига материнской функции, соответственно. Под вейвлет-функциями в рамках настоящей работы понимается хорошо известное семейство вейвлетов Добеши.
Начнем анализ полученных результатов со сравнения классической и "гистерезисной" балок с гистерезисным преобразователем Боука — Вена. Как видно из рис. 16 на начальных временных интервалах (с 0 по 3 отчет модельного времени) балки ведут себя одинаково, однако спустя 8 интервала временного отчета "гистерезисная" балка стаблизируется, в то время как классическая балка демонстрирует периодические колебания. Такая структура колебаний объясняется диссипирующими свойствами гистерезисного преобразователя (рис. 17), который трансформирует энергию колебаний сейсмической волны в иные формы.
Также из результатов моделирования, представленных на рис. 16, видно, что реация классической и "гистерезисной" балок различна. А именно в моментах воздействия сейсмической волны гистерезисная балка демонстрирует возмущения меньшой амплитуды (особенно четко это наблюдается для балки с использованием модели Прандтля — Ишлинского) в сравнении с классической балкой. После прекращения сейсмического воздействия обе балки стабилизируются в состояние равновесия, однако гистерезисная балка, формализуемая преобразователем Прандтля — Ишлинского, оказывается в положении равновесия, отличном от исходного (рис. 16, б), что согласуется с упруго-пластическими свойствами гистерезисного материала.
Также, как видно из представленных результатов, обе балки все еще продолжают совершать колебания после воздействия нагрузки. Однако балка с учетом гистерезисных свойств показывает колебания значительно меньшей амплитуды, чем классическая (0.0004 против 0.1 соответственно).
Из продемонстрированных результатов можно сделать вывод, что учет гистерезисных свойств позволяет идентифицировать новые формы поведения в сравнении с классической моделью.
х
1.0
0.5
0.0
-0.5
-1.0
10 t
15
20
б
Рис. 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
0.0
rUxx(L/2, t)
а
Рис. 17. Параметрические графики гистерезисных петель (время t как параметр): входно-выходная зависимость входа uxx от выхода Г[ихх] в точке x = L/2: а — преобразователь Прандтля — Ишлинского,
б — преобразователь Боука — Вена Fig. 17. Parametric graphs of hysteresis loops (time t as a parameter): input-output dependence of the input uxx on the output Г[ихх] at the point x = L/2: a — the Prandtl — Ishlinskii operator, b — the Bouc — Wen
operator
а
Заключение
В настоящей статье исследована динамика балки с распределенными гистерезисными свойствами. В первой части работы получены уравнения с использованием гамильтоновско-го подхода, описывающие модель движений балки, когда носитель гистерезисных свойств распределен по всей ее длине. Приведены результаты моделирования колебаний балки, находящейся под внешней нагрузкой различной природы: сейсмической (формализуемой материнской вейвлет-функцией из семейства вейвлетов Добеши), ударной и периодической. Произведено сравнение характера колебаний классической и "гистерезисной" балок. Для описания гистерезисных свойств использовались два подхода: феноменологический (преобразователь Боука — Вена) и конструктивный (преобразователь Прандтля — Ишлинского). На основе численных экспериментов установлено, что балка с распределенными гистерезисными свойствами обладает повышенной (по сравнению с классической моделью) устойчивостью по отношению к внешним
воздействиям различной природы. Иными словами, области устойчивости (в пространстве параметров этой системы) шире для гистерезисной балки. Полученные результаты могут найти применение при разработке и проектировании несущих конструкций в сейсмоопасных регионах.
Финансирование. Статья поддержана РНФ. Грант № 23-29-00696.
Информация о конфликте интересов: автор и рецензенты заявляют об отсутствии конфликта интересов.
Цитирование. Карпов Е.А. Динамика балки Эйлера — Бернулли с учетом распределенных гистерезисных свойств // Вестник Самарского университета. Естественнонаучная серия / Vestnik of Samara State University. Natural Science Series 2024. Т. 30, № 3. C. 35-62. DOI: 10.18287/2541-7525-2024-30-335-62.
© Карпов Е.А., 2024
Карпов Евгений Андреевич ([email protected]) - аспирант кафедры цифровых технологий, Воронежский государственный университет, 394018, Российская Федерация, г. Воронеж, Университетская пл., 1.
Литература
[1] Timoshenko S. History of Strength of Materials: With a Brief Account of the History of Theory of Elasticity and Theory of Structures. Dover Civil and Mechanical Engineering Series. New York: Dover Publications. 1983. ISBN 9780486611877. URL: https://archive.org/details/historyofstrengt0000timo_k8r2.
[2] Elishakoff I. Who developed the so-called Timoshenko beam theory? // Mathematics and Mechanics of Solids. 2020. Vol. 25, issue 1. P. 97-116. DOI: https://doi.org/10.1177/1081286519856931.
[3] Bauchau O.A., Craig J.I. Euler — Bernoulli Beam Theory // Solid Mechanics and Its Applications. Dordrecht: Springer Netherlands, 2009. P. 173-221. DOI: http://doi.org/10.1007/978-90-481-2516-6_5.
[4] Esen I. Dynamics of size-dependant Timoshenko micro beams subjected to moving loads // International Journal of Mechanical Sciences. 2020. Vol. 175. P. 105501. DOI: https://doi.org/10.1016/j.ijmecsci.2020.105501.
[5] Krysko A., Awrejcewicz J., Kutepov I., Krysko V. Stability of curvilinear Euler — Bernoulli beams in temperature fields // International Journal of Non-Linear Mechanics. 2017. Vol. 94. P. 207-215. DOI: https://doi.org/10.1016/j.ijnonlinmec.2016.12.004.
[6] Semenov M., Reshetova O., Borzunov S., Meleshenko P. Self-oscillations in a system with hysteresis: the small parameter approach // The European Physical Journal Special Topics. 2021. Vol. 230. P. 3565-3571. DOI: https://doi.org/10.1140/epjs/s11734-021-00237-3.
[7] Semenov M., Solovyov A., Meleshenko P., Balthazar J. Nonlinear Damping: From Viscous to Hysteretic Dampers // Recent Trends in Applied Nonlinear Mechanics and Physics. Cham: Springer, 2018. P. 259-275. DOI: https://doi.org/10.1007/978-3-319-63937-6_15.
[8] Semenov M., Solovyov A., Meleshenko P., Reshetova O. Efficiency of hysteretic damper in oscillating systems // Mathematical Modelling of Natural Phenomena. 2020. Vol. 15, article number 43. DOI: https://doi.org/10.1051/mmnp/2019053.
[9] Semenov M., Reshetova O., Tolkachev A., Solovyov A., Meleshenko P. Oscillations Under Hysteretic Conditions: From Simple Oscillator to Discrete Sine-Gordon Model // Belhaq M. (eds.) Topics in Nonlinear Mechanics and Physics. Springer Proceedings in Physics. Vol. 228. Singapore: Springer, 2019. P. 229-253. DOI: https://doi.org/10.1007/978-981-13-9463-8_12.
[10] Medvedskii A., Meleshenko P., Nesterov V., Reshetova O., Semenov M., Solovyov A. Unstable oscillating systems with hysteresis: Problems of stabilization and control // Journal of Computer and Systems Sciences International. 2020. Vol. 59. P. 533-556. DOI: http://dx.doi.org/10.31857/S0002338820030099.
[11] Semenov M., Borzunov S., Meleshenko P., Stochastic preisach operator: definition within the design approach // Nonlinear Dynamics. 2020. Vol. 101. P. 2599-2614. DOI: https://doi.org/10.1007/s11071-020-05907-w.
[12] Semenov M., Borzunov S., Meleshenko P. A new way to compute the Lyapunov characteristic exponents for non-smooth and discontinues dynamical systems // Nonlinear Dynamics. 2022. Vol. 109. P. 1805-1821. DOI: https://doi.org/10.1007/s11071-022-07492-6.
[13] Semenov M., Meleshenko P., Borzunov S., Reshetova O., Barsukov A. A Simple Model of the Energy Harvester within a Linear and Hysteresis Approach // Micromachines. 2023. Vol. 14, no. 2, article number 310. DOI: https://doi.org/10.3390/mi14020310. EDN: https://elibrary.ru/wssbtt.
[14] Борзунов С. Трансформация колебаний неустойчивой системы в преобразователе-накопителе энергии // Вестник Самарского университета. Естественнонаучная серия / Vestnik of Samara University. Natural Science Series. 2023. Т. 29, № 2. С. 7-18. DOI: http://dx.doi.org/10.18287/2541-7525-2023-29-2-7-18.
[15] Semenov M., Solovyov A., Meleshenko P. Stabilization of coupled inverted pendula: From discrete to continuous case // Journal of Vibration and Control. 2021. Vol. 27, issue 1-2. P. 43-56. DOI: https://doi.org/10.1177/1077546320923436.
[16] Osintsev M., Sobolev V. Order Reduction of Kalman-Bucy Filter for Systems with Low Measurement Noise. Slow-Fast Systems and Hysteresis: Theory and Applications // Korobeinikov A. (eds.) Extended Abstracts Summer 2016. Trends in Mathematics. Vol. 10. Birkhauser, Cham, 2018. P. 47-52. DOI: https://doi.org/10.1007/978-3-030-01153-6_9.
[17] Sobolev V. Thrice Critical Case in Singularly Perturbed Control Problems: Slow-Fast Systems and Hysteresis: Theory and Applications // Korobeinikov A. (eds.) Extended Abstracts Summer 2016. Trends in Mathematics. Vol. 10. Birkhauser, Cham, 2018. P. 83-87. DOI: https://doi.org/10.1007/978-3-030-01153-6_15.
[18] Tollmien W., Schlichting H., Gortler H., Riegels F. Ein Gedankenmodell zur kinetischen Theorie der festen Korper // Riegels F.W. (eds.) Ludwig Prandtl Gesammelte Abhandlungen. Berlin; Heidelberg: Springer, 1961. P. 149-184. DOI: https://doi.org/10.1007/978-3-662-11836-8_12.
[19] Preisach F. Uber die magnetische nachwirkung // Zeitschrift fur Physik. 1935. Vol. 94, issue 5. P. 277-302. DOI: https://doi.org/10.1007/BF01349418.
[20] Iwan W. A Distributed-Element Model for Hysteresis and Its Steady-State Dynamic Response // Journal of Applied Mechanics. 1966. Vol. 33, issue 4. P. 893-900. DOI: http://doi.org/10.1115/1.3625199.
[21] Bouc R. Forced vibrations of mechanical systems with hysteresis // Proceedings of the fourth Conference on Nonlinear Oscillations, Prague, September 5-9. 1967. P. 315-321.
[22] Bouc R. Modеle ma^matique d'hysteresis // Acustica. 1971. Vol. 24. P. 16-25.
[23] Lin Y., Cai G. Random vibration of hysteretic systems // Schiehlen W. (eds.) Nonlinear Dynamics in Engineering Systems. International Union of Theoretical and Applied Mechanics. Berlin; Heidelberg: Springer, 1976. P. 189-196. DOI: https://doi.org/10.1007/ 978-3-642-83578-0_24.
[24] Krasnoselskii M., Niezgodka M., Pokrovskii A. Systems with Hysteresis. Berlin; Heidelberg: Springer, 2012. 410 p. ISBN 978-3-642-61302-9. DOI: https://doi.org/10.1007/978-3-642-61302-9.
[25] Visintin A. Chapter 1 - mathematical models of hysteresis. In: Bertotti G., Mayergoyz I.D. (eds.) The Science of Hysteresis. Oxford: Academic Press, 2006. ISBN 978-0-12-480874-4. P. 1-123. URL: https://www.science.unitn.it/ visintin/Elsevier2006.pdf.
[26] Desch W., Turi J. The stop operator related to a convex polyhedron // Journal of Differential Equations. 1999. Vol. 157, issue 2. P. 329-347.
[27] Lang H., Dressler K., Pinnau R., Speckert M. Notes on Lipschitz estimates for the stop and play operator in plasticity // Applied Mathematics Letters. 2009. Vol. 22, issue 4. P. 623-627. DOI: https://doi.org/10.1016/j.aml.2008.07.004.
[28] Matsuo T., Terada Y., Shimasaki M. Representation of minor hysteresis loops of a silicon steel sheet using stop and play models // Physica B: Condensed Matter. 2006. Vol. 372, issues 1-2. P. 25-29. DOI: https://doi.org/10.1016/j.physb.2005.10.121.
[29] Al Janaideh M., Krejch P. An inversion formula for a Prandtl-Ishlinskii operator with time dependent thresholds // Physica B: Condensed Matter. 2011. Vol. 406, issue 8. P. 1528-1532. https://doi.org/10.1016/j.physb.2011.01.062.
[30] Li Z., Zhang X., Ma L. Development of a combined Prandtl-Ishlinskii-Preisach model // Sensors and Actuators A: Physical. 2020. Vol. 304. P. 111797. DOI: https://doi.org/10.1016/j.sna.2019.111797.
[31] Krejch P., Sprekels J. On a class of multi dimensional Prandtl-Ishlinskii operators // Physica B: Condensed Matter. 2001. Vol. 306, issues 1-4, P. 185-190. DOI: https://doi.org/10.1016/S0921-4526(01)01001-8.
[32] Hassani V., Tjahjowidodo T., Do T.N. A survey on hysteresis modeling, identification and control // Mechanical Systems and Signal Processing. 2014. Vol. 49, issues 1-2. P. 209-233. DOI: https://doi.org/10.1016/j.ymssp.2014.04.012.
[33] Dargahi A., Rakheja S., Sedaghati R. Development of a field dependent Prandtl-Ishlinskii model for magnetorheological elastomers // Materials & Design. 2019. Vol. 166. P. 107608. DOI: https://doi.org/10.1016/j.matdes.2019.107608.
[34] Krejch P. Reliable solutions to the problem of periodic oscillations of an elastoplastic beam // International Journal of Non-Linear Mechanics. 2002. Vol. 37, issue 8. P. 1337-1349. DOI: https://doi.org/10.1016/S0020-7462(02)00022-7.
[35] Ikhouane F., Rodellar J. On the hysteretic Bouc-Wen model // Nonlinear Dynamics. 2005. Vol. 42, issue 1. P. 63-78. DOI: https://doi.org/10.1007/s11071-005-0069-3.
[36] Ikhouane F., Rodellar J., On the hysteretic Bouc-Wen model // Nonlinear Dynamics. 2005. Vol. 42, issue 1. P. 79-95. DOI: https://doi.org/10.1007/s11071-005-0070-x.
[37] Ikhouane F., Hurtado J.E., Rodellar J. Variation of the Hysteresis Loop with the Bouc-Wen Model Parameters // Nonlinear Dynamics. 2007. Vol. 48, issue 4. P. 361-380. DOI: https://doi.org/10.1007/s11071-006-9091-3.
MA THEMA TICAL MODELLING
[38] Ismail M., Ikhouane F., Rodellar J. The Hysteresis Bouc-Wen model, a Survey // Archives of Computational Methods in Engineering. 2009. Vol. 16, issue 2. P. 161-188. DOI: https://doi.org/10.1007/s11831-009-9031-8.
[39] Semenov M., Karpov E., Tikhomirov S., Meleshenko P., Teplyakova M. Stabilization of Chaos Via Strong Nonlinearities: The Lorenz-Malkus Wheel Under Coulomb and Hystersis Frictions // Balthazar J.M. (eds.) Vibration Engineering and Technology of Machinery. Mechanisms and Machine Science. Vol. 95. Cham: Springer, 2021. P. 3-36. DOI: https://doi.org/10.1007/978-3-030-60694-7_1.
DOI: 10.18287/2541-7525-2024-30-3-35-62
Dynamics of the Euler — Bernoully beam with distributed hysteresis properties
Karpov E.A.
Voronezh State University, Voronezh, Russian Federation; [email protected];
Abstract. In this paper, we present a new mathematical approach to the analysis of a beam with distributed hysteresis properties. These hysteresis characteristics are described by two methods: phenomenological (Bouc — Wen model) and constructive (Prandtl — Ishlinskii model). The equations for beam are developed using the well-known Hamilton method. We investigate the dynamic response of a hysteresis beam under various external loads, including impulse, periodic and seismic loads. The results of numerical simulations show that the hysteresis beam exhibits differently to external influences as compared to the classical Euler-Bernoulli beam. In particular, under the same external loads, the vibration amplitude and energy characteristics of the hysteresis beam are lower than those of the classical one. These findings can be useful for buildings developers in the design of external load resistant buildings and structures.
Key words: hysteresis; Euler-Bernoulli beam; Bouc — Wen model; Prandtl — Ishlinskii model; nonlinear dynamics; stability; elastoplasticity.
Funding. The article is supported by the RNF. Grant № 23-29-00696.
Information about the conflict of interests: the author and reviewers declare no conflicts of interest
Citation. Karpov E.A. Dynamics of the Euler — Bernoully beam with distributed hysteresis properties. Vestnik Samarskogo universiteta. Estestvennonauchnaya seriya / Vestnik of Samara State University. Natural Science Series, 2024, vol. 30, no. 3, pp. 35-62. DOI: 10.18287/2541-7525-2024-30-3-3562. (In Russ.)
© Karpov E.A., 2024
Evgeny A. Karpov ([email protected]) - postgraduate student, student of the Department of Digital Technologies, Voronezh State University, 1, Universitetskaya Square, Voronezh, 394018, Russian Federation.
Received: 05.06.2024 Revised: 02.08.2024 Accepted: 02.09.2024
Scientific article
References
[1] Timoshenko S. History of Strength of Materials: With a Brief Account of the History of Theory of Elasticity and Theory of Structures. Dover Civil and Mechanical Engineering Series. New York: Dover Publications, 1983. ISBN 9780486611877. Available at: https://archive.org/details/historyofstrengt0000timo_k8r2.
[2] Elishakoff I. Who developed the so-called Timoshenko beam theory? Mathematics and Mechanics of Solids, 2020, vol. 25, issue 1, pp. 97-116. DOI: https://doi.org/10.1177/1081286519856931.
[3] Bauchau O.A., Craig J.I. Euler — Bernoulli Beam Theory. In: Solid Mechanics and Its Applications. Dordrecht: Springer Netherlands, 2009. ISBN 978-90-481-2516-6, pp. 173-221. DOI: http://doi.org/10.1007/978-90-481-2516-6_5.
[4] Esen I. Dynamics of size-dependant Timoshenko micro beams subjected to moving loads. International Journal of Mechanical Sciences, 2020, vol. 175, p. 105501. DOI: https://doi.org/10.1016/j.ijmecsci.2020.105501.
[5] Krysko A., Awrejcewicz J., Kutepov I., Krysko V. Stability of curvilinear Euler — Bernoulli beams in temperature fields. International Journal of Non-Linear Mechanics, 2017, vol. 94, pp. 207-215. DOI: https://doi.org/10.1016/j.ijnonlinmec.2016.12.004.
[6] Semenov M., Reshetova O., Borzunov S., Meleshenko P. Self-oscillations in a system with hysteresis: the small parameter approach. The European Physical Journal Special Topics, 2021, vol. 230, pp. 3565-3571. DOI: https://doi.org/10.1140/epjs/s11734-021-00237-3.
[7] Semenov M., Solovyov A., Meleshenko P., Balthazar J. Nonlinear Damping: From Viscous to Hysteretic Dampers. Recent Trends in Applied Nonlinear Mechanics and Physics. Cham: Springer, 2018, pp. 259-275. DOI: https://doi.org/10.1007/978-3-319-63937-6_15.
[8] Semenov M., Solovyov A., Meleshenko P., Reshetova O. Efficiency of hysteretic damper in oscillating systems. Mathematical Modelling of Natural Phenomena, 2020, vol. 15, Article number 43. DOI: https://doi.org/10.1051/mmnp/2019053.
[9] Semenov M., Reshetova O., Tolkachev A., Solovyov A., Meleshenko P. Oscillations Under Hysteretic Conditions: From Simple Oscillator to Discrete Sine-Gordon Model. In: Belhaq M. (eds.) Topics in Nonlinear Mechanics and Physics. Springer Proceedings in Physics, vol. 228. Singapore: Springer, 2019, pp. 229-253. DOI: https://doi.org/10.1007/978-981-13-9463-8_12.
[10] Medvedskii A., Meleshenko P., Nesterov V., Reshetova O., Semenov M., Solovyov A., Unstable oscillating systems with hysteresis: Problems of stabilization and control. Journal of Computer and Systems Sciences International, 2020, vol. 59, pp. 533-55б. DOI: http://dx.doi.org/10.31857/S0002338820030099.
[11] Semenov M., Borzunov S., Meleshenko P., Stochastic preisach operator: definition within the design approach. Nonlinear Dynamics, 2020, vol. 101, pp. 2599-2б14. DOI: https://doi.org/10.1007/s11071-020-05907-w.
[12] Semenov M., Borzunov S., Meleshenko P. A new way to compute the Lyapunov characteristic exponents for non-smooth and discontinues dynamical systems. Nonlinear Dynamics, 2022, vol. 109, pp. 1805-1821. DOI: https://doi.org/10.1007/s11071-022-07492-6.
[13] Semenov M., Meleshenko P., Borzunov S., Reshetova O., Barsukov A. A Simple Model of the Energy Harvester within a Linear and Hysteresis Approach. Micromachines, 2023, vol. 14, no. 2, Article number 310. DOI: https://doi.org/10.3390/mi14020310. EDN: https://elibrary.ru/wssbtt.
[14] Borzunov S. Transformation of oscillations of an unstable system in an energy harvester. Vestnik of Samara University. Natural Science Series, 2023, vol. 29, no. 2, pp. 7-18. DOI: http://dx.doi.org/10.18287/2541-7525-2023-29-2-7-18. (In Russ.)
[15] Semenov M., Solovyov A., Meleshenko P. Stabilization of coupled inverted pendula: From discrete to continuous case. Journal of Vibration and Control, 2021, vol. 27, issue 1-2, pp. 43-56. DOI: https://doi.org/10.1177/1077546320923436.
[16] Osintsev M., Sobolev V. Order Reduction of Kalman — Bucy Filter for Systems with Low Measurement Noise: Slow-Fast Systems and Hysteresis: Theory and Applications. In: Korobeinikov A. (eds.) Extended Abstracts Summer 2016. Trends in Mathematics, vol 10. Birkhauser, Cham, 2018, pp. 47-52. DOI: https://doi.org/10.1007/978-3-030-01153-6_9.
[17] Sobolev V. Thrice Critical Case in Singularly Perturbed Control Problems: Slow-Fast Systems and Hysteresis: Theory and Applications. In: Korobeinikov A. (eds.) Extended Abstracts Summer 2016. Trends in Mathematics, 2018, vol. 10. Birkhauser, Cham, pp. 83-87. DOI: https://doi.org/10.1007/978-3-030-01153-6_15.
[18] Tollmien W., Schlichting H., Görtler H., Riegels F., Ein Gedankenmodell zur kinetischen Theorie der festen Körper. In: Riegels F.W. (eds.) Ludwig Prandtl Gesammelte Abhandlungen. Berlin; Heidelberg: Springer, 1961, pp. 149-184. DOI: https://doi.org/10.1007/978-3-662-11836-8_12.
[19] Preisach F. Über die magnetische nachwirkung. Zeitschrift fur Physik, 1935, vol. 94, issue 5. pp. 277-302. DOI: https://doi.org/10.1007/BF01349418.
[20] Iwan W. A Distributed-Element Model for Hysteresis and Its Steady-State Dynamic Response. Journal of Applied Mechanics, 1966, vol. 33 (4), pp. 893-900. DOI: http://doi.org/10.1115/1.3625199.
[21] Bouc R. Forced vibrations of mechanical systems with hysteresis. Proceedings of the 4th Conference on Nonlinear Oscillations, Prague, September ,5-9, 1967, pp. 315-321.
[22] Bouc R. Modele mathematique d'hysteresis. Acustica, 1971, vol. 24, pp. 16-25.
[23] Lin Y., Cai G. Random Vibration of Hysteretic Systems. In: Schiehlen W. (eds.) Nonlinear Dynamics in Engineering Systems. International Union of Theoretical and Applied Mechanics. Berlin; Heidelberg: Springer, 1976, pp. 189-196. DOI: https://doi.org/10.1007/978-3-642-83578-0_24.
[24] Krasnoselskii M., Niezgodka M., Pokrovskii A. Systems with Hysteresis. Berlin; Heidelberg: Springer, 2012, 410 p. ISBN 978-3-642-61302-9. DOI: https://doi.org/10.1007/978-3-642-61302-9.
[25] Visintin A. Chapter 1 - Mathematical Models of Hysteresis. In: Bertotti G., Mayergoyz I.D. (eds.) The Science of Hysteresis. Oxford: Academic Press, 2006, ISBN 978-0-12-480874-4, pp. 1-123. Available at: https://www.science.unitn.it/ visintin/Elsevier2006.pdf.
[26] Desch W., Turi J. The stop operator related to a convex polyhedron. Journal of Differential Equations, 1999, vol. 157, issue 2, pp. 329-347.
[27] Lang H., Dressler K., Pinnau R., Speckert M. Notes on Lipschitz estimates for the stop and play operator in plasticity. Applied Mathematics Letters, 2009, vol. 22, issue 4, pp. 623-627. DOI: https://doi.org/10.1016/j.aml.2008.07.004.
[28] Matsuo T., Terada Y., Shimasaki M. Representation of minor hysteresis loops of a silicon steel sheet using stop and play models. Physica B: Condensed Matter, vol. 372, issues 1-2, pp. 25-29. DOI: https://doi.org/10.1016/j.physb.2005.10.121.
[29] Al Janaideh M., Krejch P. An inversion formula for a Prandtl-Ishlinskii operator with time dependent thresholds. Physica B: Condensed Matter, 2011, vol. 406, issue 8, pp. 1528-1532. DOI: https://doi.org/10.1016/j.physb.2011.01.062.
[30] Li Z., Zhang X., Ma L. Development of a combined Prandtl Ishlinskii — Preisach model. Sensors and Actuators A: Physical, 2020, vol. 304, p. 111797. DOI: https://doi.org/10.1016/j.sna.2019.111797.
[31] Krejch P., Sprekels J. On a class of multi dimensional Prandtl — Ishlinskii operators. Physica B: Condensed Matter, 2001, vol. 306, issues 1-4, pp. 185-190. DOI: https://doi.org/10.1016/S0921-4526(01)01001-8.
[32] Hassani V., Tjahjowidodo T., Do T.N. A survey on hysteresis modeling, identification and control. Mechanical Systems and Signal Processing, 2014, vol. 49, issues 1-2, pp. 209-233. DOI: https://doi.org/10.1016/j.ymssp.2014.04.012.
[33] Dargahi A., Rakheja S., Sedaghati R. Development of a field dependent Prandtl-Ishlinskii model for magnetorheological elastomers. Materials & Design, 2019, vol. 166, p. 107608. DOI: https://doi.org/10.1016/j.matdes.2019.107608.
[34] Krejch P. Reliable solutions to the problem of periodic oscillations of an elastoplastic beam. International Journal of Non-Linear Mechanics, 2002, vol. 37, issue 8, pp. 1337-1349. DOI: https://doi.org/10.1016/S0020-7462(02)00022-7.
[35] Ikhouane F., Rodellar J. On the hysteretic Bouc-Wen model. Nonlinear Dynamics, 2005, vol. 42, issue 1, pp. 63-78. DOI: https://doi.org/10.1007/s11071-005-0069-3.
[36] Ikhouane F., Rodellar J. On the hysteretic Bouc-Wen model. Nonlinear Dynamics, 2005, vol. 42, issue 1, pp. 79-95. DOI: https://doi.org/10.1007/s11071-005-0070-x.
[37] Ikhouane F., Hurtado J.E., Rodellar J. Variation of the Hysteresis Loop with the Bouc-Wen Model Parameters. Nonlinear Dynamics, 2007, vol. 48, issue 4, pp. 361-380. DOI: https://doi.org/10.1007/s11071-006-9091-3.
[38] Ismail M., Ikhouane F., Rodellar J. The Hysteresis Bouc—Wen Model, a Survey. Archives of Computational Methods in Engineering, 2009, vol. 16, issue 2, pp. 161-188. DOI: https://doi.org/10.1007/s11831-009-9031-8.
[39] Semenov M., Karpov E., Tikhomirov S., Meleshenko P., Teplyakova M. Stabilization of Chaos Via Strong Nonlinearities: The Lorenz-Malkus Wheel Under Coulomb and Hystersis Frictions.In: Balthazar J.M. (eds.) Vibration Engineering and Technology of Machinery. Mechanisms and Machine Science. Cham: Springer, 2021, vol. 95, pp. 3-36. DOI: https://doi.org/10.1007/978-3-030-60694-7 1.
