Dynamics of the Euler — Bernoully beam with distributed hysteresis properties

Abstract

Keywords

Введение
В настоящее время проектирование несущих конструкций различных инженерных соору-
жений, а также разработка адекватных физических моделей являются неотъемлемой частью
строительного процесса. Особенно остра задача выбора конструктивных составляющих зданий
стоит в городах и странах, находящихся в сейсмоопасных регионах. Одним из важнейших
этапов моделирования несущих конструкций в этом случае является идентификация форм
их поведения под воздействием внешних возбуждений (нагрузок). Обычно классификация
нагрузок выполняется по нескольким критериям: по характеру приложения (сосредоточенные и
распределенные), по продолжительности во времени (переменные и постоянные), по характеру
действия (статические и динамические). В зависимости от типа нагрузки проявляются те или
иные свойства несущих конструкций, являющиеся, в свою очередь, следствием внутренних
особенностей используемых материалов.
Timoshenko Theory Euler – Bernoulli Theory
Timoshenko Theory
Рис. 1. Графическое представление принципиального различия между теориями Тимошенко
и Эйлера — Бернулли
Fig. 1. Illustration of the main difference between Timoshenko and Euler — Bernoulli beam theories
Карпов Е.А. Динамика балки Эйлера — Бернулли с учетом распределенных гистерезисных свойств
Karpov E.A. Dynamics of the Euler — Bernoully beam with distributed hysteresis properties 36 из 62
Известно, что одним из базовых элементов любой несущей конструкции является балка.
Математическое описание балки обычно формализуется в рамках теории Эйлера — Бернулли [1]
или теории Тимошенко [2] (рис. 1). Основное различие между вышеозначенными подходами
заключается в поведении балки при поперечных деформациях. В теории Эйлера — Бернулли
предполагается, что во время деформаций поперечное сечение перпендикулярно нейтральной
оси. В то же время в теории Тимошенко оно не перпендикулярно нейтральной оси и напрямую
зависит от интенсивности деформации. В настоящей статье модель колебаний балки строится
с использованием теории Эйлера — Бернули в предположении, что балка имеет шарнирное
закрепление в двух концах и обладает нулевыми начальными условиями (балка покоится и не
имеет начальной скорости). Следуя работе [3], классическое уравнение колебаний балки имеет
следующий вид: 


utt + a2uxxxx = 1(x, t), t > 0, x ∈ [0, L],
u(x, 0) = 0, ut(x, 0) = 0,
u(0, t) = u(L, t) = S(t),
uxx(0, t) = uxx(L, t) = 0.
(1)
где u = u(x, t) – отклонение балки от положения равновесия, индексы обозначают частные
производные1, 1(x, t) есть функция, характеризующая вертикальную нагрузку (управление),
S(t) – функция, характеризующая внешнее воздействие, прикладываемое к закрепленным
концам балки, a2 = EI – произведение модуля упругости и момента инерции, t > 0 есть время
моделирования, L – длина балки.
Балки в зависимости от своих конструктивных особенностей находят свое применение в
различных элементах строительных (и не только) конструкций. Например, двутавровые балки
используются для обеспечения устойчивости перекрытий между этажами и в мостовых соору-
жениях, швеллерные же балки используются для армирования железобетонных конструкций,
а также в качестве рам, каркасов и т. д. Кроме того, балки по-разному реагируют на внешнюю
нагрузку, которая, в свою очередь, может быть импульсной, вибрационной, подвижной [4],
ударной, ветровой. При этом основная задача, решаемая при выборе типа балки применительно
к задачам строительства, заключается в обеспечении устойчивости всего сооружения по отно-
шению к потенциальным внешним воздействиям, к числу которых, в первую очередь, относятся
сейсмические волны как природного, так и техногенного характера. В таком случае балка
выступает в роли демпфирующего элемента, трансформируя входящую энергию воздействия в
иные ее формы [5].
Известно, что во время нагрузок балка начинает деформироваться, порождая во внут-
ренней структуре материала силы, стремящиеся вернуть ей исходную форму. Если после
снятия нагрузки она принимает форму, отличную от исходной, то следует говорить об упруго-
пластическом гистерезисе. Он играет важную роль в современных исследованиях реальных
систем и процессов. Это явление было обнаружено во многих областях современной науки:
физике [6–14], экономике, биологии, химии и т. д. Недавние исследования показывают важность
и целесообразность использования гистерезиса в нейронных сетях, например, включая его
как звено или как самостоятельный элемент в функцию активации соответствующей сети.
Такой широкий спектр применения обуславливается важными структурными особенностями
гистерезиса – зависимостью текущего состояния системы от предыстории и начального состоя-
ния как от параметра. Еще более важной особенностью гистерезиса является тот факт, что
в ряде задач он позволяет моделировать процесс диссипации энергии, что, в свою очередь,
позволяет использовать гистерезисные преобразователи в качестве управляющего элемента,
стабилизирующего поведение динамических систем [15–17].
Теоретические модели, описывающие гистерезисные явления, восходят к классическим
работам Прандтля [18], Прейзаха [19] и др. В начале 60-х годов XX века появились первые
феноменологические модели [20–23]. Теория конструктивных моделей гистерезиса, учитываю-
1Здесь uxxxx =
∂4u
∂x4 , utt =
∂2u
∂t2 .
Вестник Самарского университета. Естественнонаучная серия 2024. Том 3, № 30. С. 35–62
Vestnik of Samara State University. Natural Science Series 2024, vol. 3, no. 30, pp. 35–62 37 из 62
щая физические особенности конструктивных элементов исследуемых систем, получила свое
развитие в трудах М.А. Красносельского и его учеников в Воронежском университете в 60-х
годах прошлого века. Фундаментальная монография [24] с подробным описанием свойств
конструктивных моделей гистерезиса вышла в 1983 г. Разработанная в ней методология позво-
лила описывать гистерезис на языке нелинейных операторов, зависящих от своего начального
состояния как от параметра. Математический аппарат, предложенный в ней, базируется на
выделении элементарных носителей гистерезисных свойств – гистеронов – преобразователей,
наделенных пространством состояний и входно-выходными соответствиями. Такой подход
позволил имплементировать модели гистерезиса в состав более сложных систем, описываемых
дифференциальными уравнениями.
В настоящей статье предлагается модель балки с учетом распределенных гистерезисных
свойств, в разд. 1 описываются используемые в работе модели гистерезиса, разд. 2 посвящен
выводу уравнений колебаний гистерезисной балки, в разд. 3 описаны численные методы,
используемые в работе, а в разд. 4 приводятся результаты моделирования и сравнительный
анализ классической и гистерезисной балок. В качестве моделируемой нагрузки будут выступать
периодические и импульсные воздействия, а также нагрузка сейсмического характера.
1. Модели гистерезиса: феноменологический и конструктивный
подходы
1.1. Преобразователь Прандтля
Одной из самых известных и простых гистерезисных моделей является преобразователь
Прандтля [25–27] (также известный как “stop operator”). Он используется для моделирования
упруго-пластических волокон [28], где состояния полностью определяются величинами пере-
менной деформации и переменного напряжения, а также применяется в качестве элементарной
составляющей в сложных континуальных гистерезисных моделях в рамках теории упруго-
пластического гистерезиса. Действительно, преобразователь Прандтля идеально описывает
поведение материала: материал, находящийся в изначальном состоянии, когда нормальное на-
пряжение и деформация отсутствуют, под нагрузкой начинает вести себя упруго до некоторого
порогового значения h, а после демонстрирует пластические свойства, находясь под постоянным
напряжением. Это поведение также называется “пластическое течение”, так как при уменьше-
нии напряжения материал обладает пластическими свойствами, пока не достигнет порогового
значения с противоположным знаком −h. Однако отметим, что преобразователь Прандтля
позволяет моделировать упруго-пластические свойства без деформационного упрочнения2.
Структурно этот преобразователь может быть представлен как последовательное соединение
линейного упругого элемента E (например, пружина) с жестким, идеально пластическим
элементом P (например, брусок, двигающийся по поверхности, под действием силы сухого
трения) (рис. 2). Если предположить, что упругость пружины E равна единице, следует говорить
о более простом преобразователе – упоре.
Описание преобразователя Прандтля начнем с определения его состояний. Состояниями
этого преобразователя являются пары “вход-выход” {x, u} с областью определения в виде полосы
Ω = Ω(ΓP, E), угол наклона которой зависит от коэффициента E, определяющего модуль
упругости. Причем полоса Ω зависит также и от порогового значения h (полоса пересекает
ось Ox в точках ±h), которое является еще одной характеристикой этой модели. Иначе говоря,
модель может быть формализована в виде равенства:
σ(t) = ΓP[h]ε(t),
2Также известное как “strain hardening”, которое представляет собой упрочнение материала (например,
металла или полимера) за счет пластической деформации, которое появляется из-за возникновения и движения
дислокаций в кристаллической структуре материала
Карпов Е.А. Динамика балки Эйлера — Бернулли с учетом распределенных гистерезисных свойств
Karpov E.A. Dynamics of the Euler — Bernoully beam with distributed hysteresis properties 38 из 62
E P
Рис. 2. Графическое представление гистерезисного преобразователя Прандтля. Здесь E и P — упругий
и пластический элементы соответственно
Fig. 2. Graphical representation of the Prandtl operator. E and P are the elastic and the plastic elements
respectively
где
ΓP[h]ε(t) =

min{h,φ(t, τ)}, если ε(t) ≥ ε(t − τ),
max{−h,φ(t, τ)}, иначе
φ(t, τ) = E [ε(t) − ε(t − τ)] + σ(t − τ)
или в более компактной форме:
σ(t) = ΓP[h]ε(t) = min{h, max{−h,φ(t, τ)}}, (2)
где σ(t − τ) — значение выхода гистерезисного преобразователя в момент времени t − τ, ε(t) —
входной сигнал, σ(t) — выходной сигнал (реакция преобразователя). Графическое представление
описанного выше преобразователя отображено на рис. 3.
Рис. 3. График зависимости входа ε(t) от выхода σ(t) для преобразователя Прандтля с пороговым
значением, равным h, и областью определения Ω = Ω(ΓP, 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 Ω = Ω(ΓP, E) presented as a band with the slope angle determined by the coefficient E
1.2. Конструктивный подход: преобразователь Прандтля — Ишлинского
Одним из самых известных преобразователей в рамках конструктивного подхода к моде-
лированию упруго-пластического гистерезиса является преобразователь Прандтля — Ишлин-
ского. Обычно этот преобразователь применяется в задачах, где требуется формализовать
гистерезисную связь между напряжением и деформацией. Он является компиляцией преоб-
разователей Прандтля и материала Ишлинского. В отличие от преобразователя Прандтля
преобразователь Прандтля — Ишлинского позволяет учитывать деформационное упрочнение
совместно с упруго-пластическими свойствами (рис. 4). Натурные эксперименты, проведенные
с упруго-пластическими материалами, демонстрируют, что состояние материала не только
характеризуется входно-выходной характеристикой “напряжение-деформация”, но и зависит
Вестник Самарского университета. Естественнонаучная серия 2024. Том 3, № 30. С. 35–62
Vestnik of Samara State University. Natural Science Series 2024, vol. 3, no. 30, pp. 35–62 39 из 62
от некоторых внутренних переменных. Преобразователь Прандтля — Ишлинского позволяет
учитывать предысторию и является удобным инструментом для моделирования гистерезисных
зависимостей. Преобразователи с такой структурой относят к классу сложных гистерезисных
систем, описываемых посредством континуальных аналогов блок-схем (рис. 5). В настоящем
разделе приводится математическое описание этого преобразователя.
Γ
P
.......
h
1
1
Γ
P
3
Γ
P
2
h
2
h
3
h
∞
Γ
P
∞
Рис. 4. Внутренняя структура континуального преобразователя Прандтля — Ишлинского
с пороговыми значениями h = h1, h2, ..., h∞, распределенными с помощью некоторой функции
распределения и параллельно соединенными преобразователями Прандтля ΓP = ΓP1, ΓP2, ..., ΓP∞
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
ΓP = ΓP1, ΓP2, ..., ΓP∞
Зададим семейство однопараметрических преобразователей Прандтля (2), зависящих от
порогового значения h как от параметра. Пусть параметр h, характеризующий пороговое
значение, принадлежит некоторому множеству Φ. Пусть задана неубывающая непрерывная
слева функция Ξ = Ξ(h), h > 0 на множестве Φ, удовлетворяющая следующим двум условиям:
lim
h→∞
Ξ(h) = 0,
Z ∞
0
|Ξ(h)|dh < ∞.
Также пусть задано множество Z непрерывных функций σ(h), удовлетворяющих условию
|σ(h)| ≤ h.
Теперь определим множество Ω(ΓP, E) возможных состояний преобразователя ΓP как
Ω(ΓP, E) = {{ε, σ1(h)}, {ε, σ2(h)}, ..., {ε, σ∞(h)}},
где ε — это вход, σi(h) — выход (реакция) i-го преобразователя ΓP
i (каждый преобразователь
зависит от порогового значения как от параметра). Входно-выходные соответствия такого
преобразователя описываются следующим равенством:
σ(t) = ΓPI[t0, σ0(h); Ξ]ε(t) =
Z ∞
0
ΓP[t0, σ0(h); h]ε(t)dΞ(h), (3)
а переменное состояние {ε(t), σ(h, t)} определяется равенством
σ(h; t) = Q[t0, z0(h); Ξ]ε(t) = ΓP[t0, σ0(h); h]ε(t), (4)
где Q есть пучок гистеронов (конечный или в общем случае бесконечный набор гистеронов),
преобразователь ΓP описывает параллельное соединение бесконечного набора преобразователей,
Карпов Е.А. Динамика балки Эйлера — Бернулли с учетом распределенных гистерезисных свойств
Karpov E.A. Dynamics of the Euler — Bernoully beam with distributed hysteresis properties 40 из 62
описываемых соотношением (2) (см. рис. 4), σ0(h) стоит понимать как значение выхода в
предыдущий момент времени (альтернативная запись σ(t − τ, h)).
Согласно определениям и свойствам, изложенным в работе [24], подобный преобразователь
обладает свойствами статичности и виброкорректности, однако свойствами управляемости и
детерминированности не располагает. В ряде задач, вместо континуального преобразователя
Рис. 5. Блок-схема параллельного соединения упоров в модели Прандтля — Ишлинского
Fig. 5. Schema of parallel connection of the stop-operator models in the Prandtl — Ishlinskii model
Прандтля — Ишлинского удобно использовать его конечномерный аналог (именно такая ситуа-
ция имеет место в компьютерном моделировании гистерезисных зависимостей). Более того, в
моделях реальных физических систем практически всегда реализуется конечномерный аналог
континуального преобразователя. Пусть задан конечный набор размера N ∈ Z гистеронов ΓP,
формализуемых посредством преобразователя Прандтля. Предполагая, что гистероны соеди-
нены параллельно (как и в континуальной модели), и переходя от интеграла к интегральной
сумме, преобразователь примет вид (рис. 5):
σ(t) =
XN
i=1
f (hi)ΓP[t0, σi; hi]ε(t). (5)
Большое разнообразие работ, использующих преобразователь Прандтля — Ишлинского,
можно условно разделить на две категории: фундаментальные и прикладные. В рамках фунда-
ментальных работ исследуются структурные особенности преобразователя [29–31], вопросы
корректности, а также проблемы идентификации его параметров [32]. В работах приклад-
ной направленности преобразователь используется для описания гистерезисных явлений в
механических системах, например актуаторы, магнитореологические эластомеры [33], несущие
конструкции [34] и т. д.
1.3. Феноменологический подход: преобразователь Боука — Вена
Наряду с конструктивным, не менее известным является феноменологический подход к
моделированию гистерезиса. Он позволяет моделировать петли гистерезиса разной структуры,
не принимая во внимание физическую основу моделируемого процесса. Обычно феноменоло-
гический подход применяется в том случае, когда необходимо оценить реакцию системы на
сигнал лишь на основе экспериментальных данных о характеристиках входного воздействия
и отклике системы на него. Одним из самых известных преобразователей в рамках этого
подхода является преобразователь Боука — Вена, предложенный в [21; 22]. Входные-выходные
соответствия преобразователя Боука — Вена формализуются уравнением:
˙Γ
BW = x˙

A − [βsign(x˙ΓBW) + γ] |ΓBW|n    
, (6)
Вестник Самарского университета. Естественнонаучная серия 2024. Том 3, № 30. С. 35–62
Vestnik of Samara State University. Natural Science Series 2024, vol. 3, no. 30, pp. 35–62 41 из 62
где x = x(t) — входной сигнал, ΓBW = ΓBW(t) — выходной сигнал (реакция преобразователя),
γ, n, β, A > 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) =
Z t1
t0
Ldt, (7)
где t0 ≤ t ≤ t1 — некоторый промежуток времени, u = u(x, t) — кривая, описывающая динамику
Карпов Е.А. Динамика балки Эйлера — Бернулли с учетом распределенных гистерезисных свойств
Karpov E.A. Dynamics of the Euler — Bernoully beam with distributed hysteresis properties 42 из 62
балки, L — функция Лагранжа, равная разности кинетической и потенциальной энергии:
L = T −U, (8)
где U — потенциальная энергия, T — кинетическая энергия. Определим потенциальную энергию
балки как
U =
1
2
Z l
0
Mχdx, (9)
где M — изгибающий момент, а χ — деформация балки в точке x. Деформацию определим
через кривизну.
Пусть кривая, описывающая положение балки в плоскости Oxy, задана уравнением F =
= F(x, y) = 0 (F ∈ R2), тогда ее кривизна будет описываться уравнением:
κ =
|F2
yFxx − FxFyFxy + F2x
Fyy|
(F2x
+ F2
y)3/2 .
Так как продольными колебаниями можно пренебречь (согласно предположению), то урав-
нение, описывающее положение балки, будет иметь вид F(x, y) = y − u(x) = 0, а ее кривизна
определяться следующим соотношением:
κ =
uxx
(1 − u2x
)3/2 .
Ограничиваясь рассмотрением лишь малых углов отклонения поперечного сечения ϕ ≪ 1
(ux ≪ 1 и как следствие1 − u2x
−→
ux→0
1), кривизна будет определяться равенством:
κ = uxx.
В свою очередь, изменение кривизны балки равно3
χ = k − k0 = k − 0 = k = uxx (10)
где k0 — начальная кривизна покоящейся балки, равная нулю.
Для вывода уравнения изгибающего момента выберем два параллельных сечения Ω(x) и
Ω(x + Δx) на расстоянии Δx. Очевидно, что углы поворота в этих сечениях будут различными
во время деформации. Вследствие этого происходят растяжение (в верхней части балки) и
сжатие (в нижней части балки) материала балки в продольном направлении. Рассмотрим часть
балки, находящейся на расстоянии r = r(y) от ее оси симметрии. Как видно на рис. 7, длина
растянутого материала изменилась на Δs = rΔϕ. Относительное удлинение материала
ε(x, Δx, y, t) = −
Δs
Δx
= −r
Δϕ
Δx
, (11)
а деформация растяжения-сжатия в точке x
ε(x, y, t) = − lim
Δx→0
ε(x, Δx, y, t) = −r lim
Δx→0
Δϕ
Δx
= −rϕx. (12)
Угол отклонения поперечного сечения ϕ от “нормального” положения в точке x (см. рис. 6) с
учетом ортогональности оси симметрии направлению деформации равен:
ϕ = ux. (13)
3Кривизна равна нулю только в том случае, когда каждая точка балки неподвижна, однако как только балка
подвергается деформациям, то есть некоторые бесконечно малые элементы балки начинают отклоняться от
положения равновесия, кривизна становится ненулевой.
Вестник Самарского университета. Естественнонаучная серия 2024. Том 3, № 30. С. 35–62
Vestnik of Samara State University. Natural Science Series 2024, vol. 3, no. 30, pp. 35–62 43 из 62
Рис. 7. Деформированная часть балки, полученная двумя параллельными сечениями Ω(x)
и Ω(x + Δx)
Fig. 7. Deformed part of the beam obtained by means of two parallel cross sections Ω(x) and Ω(x + Δx)
Таким образом, подставляя (13) в (14), получим
ε(x, y, t) = −ruxx. (14)
Изгибающий момент, возникающий относительно оси Oz вследствие нормального напряжения
в сечении Ω(x), описывается уравнением
M(x, t) =
Z
Ω
rσdS. (15)
Теперь вычислим нормальное напряжение в сечении Ω(x), полагая, что материал балки обладает
упруго-пластическими свойствами с деформационным упрочнением. Для этого воспользуемся
преобразователем Прандтля — Ишлинского:
σ(t) = ΓPI[t0, σ0; Ξ]ε(t) :=
Z ∞
0
ΓP[t0, σ0; h]ε(t)dΞ(h), (16)
где σ(t) — выход (реакция преобразователя), ΓPI — преобразователь Прандтля — Ишлинского,
ε(t) — вход. Выполним переход от одномерного преобразователя Прандтля — Ишлинского
к двумерному, введя пространственный параметр ζ, принадлежащий множеству измеримых
параметров, по Лебегу, как предлагается в работе [34]. Тогда параметрический преобразователь
примет вид:
σ(ζ, t) = ΓPI[t0, σ0; Ξ]ε(ζ, t) :=
Z ∞
0
ΓP[t0, σ0; h]ε(ζ, t)dΞ(ζ, h),
где ζ — параметр, характеризующий пространственную (по материалу балки) плотность
распределения элементарных носителей гистерезиса. Или, используя следующее соотношение
ζ = (x, y) (параметр ζ представлен в виде пространственной переменной), окончательно
получим:
σ(x, y, t) = ΓPI[t0, σ0; Ξ]ε(x, y, t) :=
Z ∞
0
ΓP[t0, σ0; h]ε(x, y, t)dΞ(x, y, h). (17)
С учетом соотношения (17) уравнение (15) примет вид:
M(x, t) =
Z
Ω
r
 Z ∞
0
ΓP[t0, σ0; h]ε(x, y, t)dΞ(x, y, h)
!
dS. (18)
Карпов Е.А. Динамика балки Эйлера — Бернулли с учетом распределенных гистерезисных свойств
Karpov E.A. Dynamics of the Euler — Bernoully beam with distributed hysteresis properties 44 из 62
Раскрыв скобки и учитывая соотношение (14), получим
M(x, t) =
Z
Ω
Z ∞
0
rΓP[t0, σ0; h](−ruxx)dΞ(x, y, h)dS. (19)
Теперь воспользуемся свойством преобразователя Прандтля [34], а именно соотношением
rΓP[h](−ruxx) = −r2Γ[h/|r|]uxx. (20)
Для подынтегрального выражения из соотношения (19) с учетом соотношения (14)
M(x, t) =
Z
Ω
−r2
 Z ∞
0
ΓP[t0, σ0; h/|r|]uxxdΞ(x, y, h)
!
dS (21)
и введя переменнуюbΞ (x, h) =
R
Ω r2Ξ(x, y, |r|h)dS, получим
M(x, t) = −
Z ∞
0
ΓP[t0, σ0; h]uxxdhbΞ (x, h). (22)
Как можно заметить, уравнение (22) совпадает с уравнением (16), что позволяет записать
соотношение в следующей форме:
M(x, t) = −ΓPIuxx. (23)
Таким образом, связывая отношение деформации к напряжению посредством преобразователя
Прандтля — Ишлинского (17), по сути, осуществляется распределение гистерезисных свойств
вдоль длины балки.
Окончательно принимая во внимание соотношения (10) и (23), уравнение для потенциальной
энергии примет вид
U = −
1
2
Z l
0
uxxΓPI[uxx]dx. (24)
Кинетическую энергию определим классическим образом:
T =
1
2
Z l
0
ρ(x)Su2t
dx. (25)
Тогда интеграл действия для такой системы примет вид
I(u) =
Z t2
t1
1
2
 Z L
0
uxxΓPI[uxx]dx +
Z L
0
ρSu2t
dx
!
dt. (26)
Из принципа наименьшего действия в форме Гамильтона следует, что функция u = u(x, t)
является экстремалью функционала (26). Рассмотрим вариацию этой функции bu(x, t) = u(x, t)+
+ϵη(x, t), для которой справедливо следующее:
• в начальный и конечный моменты времени функция u(x, t) не варьируется;
• η(x, t1) = η(x, t2) = 0 для таких x, для которых верно 0 ≤ x ≤ L.
Таким образом, функционал (26) для кривой bu примет вид
I =
Z t2
t1
1
2
 Z L
0
(uxx + ϵηxx)Γ[uxx + ϵηxx]dx + ρS
Z L
0
(ut + ϵηt)2dx
!
dt. (27)
Ясно, что функционал действия (27) достигает экстремума в точке ϵ = 0. Согласно необ-
ходимому условию существования экстремума, положим вариацию функционала действия
относительно ϵ равной нулю
dI
dϵ
     
ϵ=0
= 0. (28)
Вестник Самарского университета. Естественнонаучная серия 2024. Том 3, № 30. С. 35–62
Vestnik of Samara State University. Natural Science Series 2024, vol. 3, no. 30, pp. 35–62 45 из 62
Перед вычислением производной по ϵ раскроем скобки внутри подынтегральных уравнений (27):
Z L
0
uxxΓPI[uxx + ϵηxx] + ϵηxxΓPI[uxx + ϵηxx]dx, (29)
ρS
Z L
0
(u2t
+ 2ϵutηt + η2t
)dx. (30)
Подставляя в (28) выражение интеграла действия (27) с учетом (29) и (30), вычислим произ-
водную по ϵ: Z t2
t1
((I1 + I2) + ρSI3)ϵ dt = 0, (31)
где
I1 =
Z L
0
uxxΓPI[uxx + ϵηxx]dx, (32)
I2 =
Z L
0
ϵηxxΓPI[uxx + ϵηxx]dx, (33)
I3 =
Z L
0
(u2t
+ 2ϵutηt + η2t
)dx. (34)
Предположим, что функция η(x, t) на концах равняется нулю, то есть η(0, t) = η(L, t) = 0.
Рассмотрим более детально каждое слагаемое (32)–(34).
Найдем частную производную по ϵ для (32). Для этого продифференцируем подынтеграль-
ное слагаемое4: Z L
0
uxx (ΓPI[uxx + ϵηxx])ϵ dx. (35)
Для дальнейших вычислений предположим, что весовая функцияbΞ дифференцируема и ее
производная есть функция f (h). Следуя этому, преобразователь Прандтля — Ишлинского
примет вид:
ΓPI[h]ϵ =
Z ∞
0
f (h)ΓP[h]ϵdh. (36)
Учитывая (36), запишем полную форму производной преобразователя Прандтля — Ишлинского
из соотношения (35):
(ΓPI[uxx + ϵηxx])ϵ =
Z ∞
0
f (h) (ΓP[uxx + ϵηxx])ϵ dh, (37)
где ΓP — преобразователь Прандтля. Представляя преобразователь Прандтля в форме (2):
ΓP[ϕ(x, t, ϵ)] = min{h, max{−h,ψ(x, t, t0, ϵ)}}, (38)
где ψ(x, t, t0, ϵ) = ΓP0 + ϕ(x, t, ϵ) − ϕ(x, t0, ϵ) и ϕ(x, t, ϵ) = uxx + ϵηxx, его производная примет
вид:
ΓP[ϕ(x, t, ϵ)]ϵ =

0, если ψ(x, t, t0, ϵ) ∈ [−h, h],
ηxx(x, t) − ηxx(x, t0), иначе.
(39)
Полученный выше результат можно представить в виде
ΓP[ϕ(x, t, ϵ)]ϵ = (ηxx(x, t) − ηxx(x, t0))×
× (Θ(ψ(x, t, t0, ϵ) + h) − Θ(ψ(x, t, t0, ϵ) − h)), (40)
4Очевидно, что производную необходимо взять только у слагаемого ΓPI[uxx + ϵηxx], поскольку u не зависит
от ϵ.
Карпов Е.А. Динамика балки Эйлера — Бернулли с учетом распределенных гистерезисных свойств
Karpov E.A. Dynamics of the Euler — Bernoully beam with distributed hysteresis properties 46 из 62
где Θ(·) — функция Хевисайда. Тогда, подставляя полученный результат в (37) и вынося
слагаемые, не зависящие от h, получим
(ΓPI[uxx + ϵηxx])ϵ = (ηxx(x, t) − ηxx(x, t0))×
×
Z ∞
0
f (h)(Θ(ψ(x, t, t0, ϵ) + h) − Θ(ψ(x, t, t0, ϵ) − h))dh. (41)
Подставив в (35) выражение (41), обозначив интеграл как некоторую функцию β(x, t, t0, ϵ):
Z L
0
uxx(ηxx(x, t) − ηxx(x, t0))β(x, t, t0, ϵ)dx, (42)
проинтегрируем два раза по частям, принимая за v функцию от η (чтобы избавиться от двух
производных по x), и окончательно получим
(I1)ϵ =
Z L
0
(η(x, t) − η(x, t0))α(x, t, t0, ϵ)dx, (43)
где α(x, t, t0, ϵ) — некоторая произвольная функция, полученная в ходе интегрирования.
Теперь рассмотрим второй интеграл I2. Дифференцируя интеграл по ϵ:
Z L
0
ηxxΓPI[uxx + ϵηxx]dx +
Z L
0
ϵηxx(ΓPI[uxx + ϵηxx])ϵdx (44)
и полагая ε = 0, получим следующее выражение:
Z L
0
ηxxΓPI[uxx]dx. (45)
Теперь продифференцируем выражение два раза по частям и окончательно получим
(I2)ϵ =
Z L
0
ηΓPI[uxx]xxdx. (46)
Нетрудно заметить, что в интеграле (34) только одно слагаемое зависит от ε. Следовательно,
дифференцируя по ε, имеем: Z L
0
2utηtdx. (47)
Теперь, продифференцировав один раз по частям полученное выражение, окончательно имеем:
(I3)ϵ =
Z L
0
2ηuttdx. (48)
Таким образом, подставляя в (31) выражения (43), (46), (48), получим:
Z t2
t1
 Z L
0
(η(x, t) − η(x, t0))α(x, t, t0, 0)dx +
Z L
0
ηΓ[uxx]xxdx
!
+
+ ρSuttηdxdt = 0. (49)
Заметим, что вышеполученное уравнение есть необходимое условие Эйлера локального экстре-
мума функционала действия. Воспользуемся леммой Дебуа-Реймона, которая является основой
для построения теории обобщенных функций, применительно к уравнению (49) и получим
окончательно уравнение колебаний гистерезисной балки5:
Γ[uxx]xx + ρSutt = 0. (50)
5Вследствие использования леммы первое интегральное слагаемое будет равно нулю их-за разности
η(x, t) − η(x, t0).
Вестник Самарского университета. Естественнонаучная серия 2024. Том 3, № 30. С. 35–62
Vestnik of Samara State University. Natural Science Series 2024, vol. 3, no. 30, pp. 35–62 47 из 62
Учитывая потенциал внешней нагрузки и полагая произведение распределением массы на
площадь поперечного сечения, равного единице, перепишем уравнение в виде
utt + Γ[uxx]xx = 1(x, t), (51)
где 1(x, t) — функция внешней нагрузки.
3. Численные методы: явная разностная схема
и метод Рунге — Кутты
Задача нахождения решения дифференциального уравнения (в общем случае системы
дифференциальных уравнений) является самой первой задачей в рамках исследования любой
системы, описываемой дифференциальными уравнениями. Как правило, обычно использу-
ются точные аналитические методы решения. Ярким примером служит применение метода
разделения переменных Фурье в задаче о колебаниях струны. Преимущество точных методов
состоит в том, что они позволяют получить решение уравнения в виде некоторой комбинации
элементарных функций (иногда решение представляется в виде квадратур от элементарных
функций). Однако большинство практических задач решить с использованием этих методов не
является возможным. Вызвано это либо особенностями области, в которой ищется решение
уравнения (многие практические задачи накладывают ограничения на пространство парамет-
ров и области решений), либо вообще невозможностью нахождения аналитического решения.
В таком случае удобно использовать численные методы, позволяющие искать решение, ап-
проксимируя исходное уравнение. Математический аппарат численных методов достаточно
развит для его применения во многих задачах. Однако применение этих методов становится
нетривиальной задачей, когда уравнения содержит операторные нелинейности гистерезисного
типа. Для решения таких уравнений необходимо модифицировать классические методы.
В рамках настоящей статьи рассматриваются численные методы, преобразующие диффе-
ренциальное уравнение в частных производных с операторной гистерезисной нелинейностью
к конечной системе алгебраических уравнений. Наиболее часто в литературе используется
следующая интерпретация разностных схем: явные и неявные разностные схемы. Неявная
разностная схема предполагает преобразование исходного уравнения к системе алгебраиче-
ских уравнений с трехдиагональной матрицей с последующим решением (например, методом
прогонки). Явная схема более проста в реализации и позволяет вычислять значение функции,
опираясь на данные, полученные на предыдущем шаге. Неявные схемы, как правило, более
устойчивы, чем явные. Однако, полагая отношение шага по времени много больше шага по
пространству, можно добиться необходимой устойчивости. Ввиду сложности формализации
гистерезисных операторов и их множественной суперпозиции в неявной разностной схеме, в
настоящей работе используется явная разностная схема.
Введем в области Ω = {0 ≤ x ≤ L, 0 ≤ t ≤ T} равномерную сетку с шагом по пространству hx
и с шагом по времени ht, равную
hx =
L
N
, ht =
T
M
, (52)
где M — количество узлов сетки по временной координате, T — время моделирования, N —
количество узлов сетки по пространственной координате (рис. 8). Причем шаги подчиняются
соотношению ht ≪ hx для устойчивости описываемой схемы. Первые два временных шага t = 0,
t = 1 вычисляются на основе начальных условий. Используя граничные условия, нетрудно
идентифицировать все значения в узлах с координатами (xi, tj), где i = 0, 1,N−1,N−2, j = 0,M.
Для вычисления остальных значений аппроксимируем уравнение колебаний гистерезисной
балки.
Уравнения, аппроксимирующие вторые производные по времени и пространству, примут вид:
utt =
uj−1
i
− 2uj
i + uj+1
i
h2t
, (53)
Карпов Е.А. Динамика балки Эйлера — Бернулли с учетом распределенных гистерезисных свойств
Karpov E.A. Dynamics of the Euler — Bernoully beam with distributed hysteresis properties 48 из 62
x
t
0
0
1
2
M
M-1
M-2
1 2
...
...
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
uxx =
uj
i−1
− 2uj
i + uj
i+1
h2x
, (54)
где uj
i соответствует значению в области Ω с координатами (i, j). Аппроксимация нулевых
начальных и граничных условий очевидна, поэтому не приводится в настоящей статье. Аппрок-
симация функции нагрузки 1(x, t) выполнена следующим образом:
1(x, t) = 1( jhx, iht), (55)
а аппроксимация гистерезисного оператора имеет вид:
Γ[uxx]xx =
Γ[uxx]
j−1
i
− 2Γ[uxx]
j
i + Γ[uxx]
j+1
i
h2x
, (56)
и, вводя новую переменную euj
i =
uj
i−1
−2uj
i+uj
i+1
h2x, согласно (54), окночательно получим соотношение:
Γ[uxx] =
Γ[euj
i ]
j−1
i
− 2Γ[euj
i ]
j
i + Γ[euj
i ]
j+1
i
h2x
, (57)
где Γ — гистерезисный оператор. Подставляя полученные аппроксимации в уравнение (51),
получим следующее соотношение:
uj−1
i
− 2uj
i + uj+1
i
h2t
+
Γ[euj
i ]
j−1
i
− 2Γ[euj
i ]
j
i + Γ[euj
i ]
j+1
i
h2x
= 1( jhx, iht). (58)
Оставляя в левой части уравнения только неизвестные, окончательно получим уравнение
колебаний гистерезисной балки, аппроксимированное явной численной схемой:
uj+1
i = 2uj
i + −uj−1
i h2t
1( jhx, iht) − h2t
Γ[euj
i ]
j−1
i
− 2Γ[euj
i ]
j
i + Γ[euj
i ]
j+1
i
h2x
. (59)
Таким образом, значение на j+1 временном слое будет высчитываться на основе j и j − 1 слоев.
Вестник Самарского университета. Естественнонаучная серия 2024. Том 3, № 30. С. 35–62
Vestnik of Samara State University. Natural Science Series 2024, vol. 3, no. 30, pp. 35–62 49 из 62
В зависимости от формализации гистерезисного звена будут рассматриваться различные
подходы к его разностностной дискретизации. Так, для расчета реакции преобразователя
Прандтля — Ишлинского будет использоваться его дискретный аналог (5). Выход преобра-
зователя Боука — Вена рассчитывается, используя модификацию хорошо известного метода
Рунге — Кутты четвертого порядка. Кратко опишем используемый алгоритм.
Пусть F = F(ξ, η) есть правая часть дифференциального уравнения преобразователя Боука —
Вена со входом ξ и значением этой функции в предыдущий момент времени η. Введем вектор
Fj =

F(ξj
2, ηj
2)
F(ξj
3, ηj
3)
...
F(ξj
N−2, ηj
N−2)

(60)
преобразователей Боука — Вена со значениями ξj
i =
uj
i−1
−2uj
i+uj
i+1
h2x
и ηj
i = Γj+1
i+1 соответственно.
Тогда модель Рунге — Кутты 4-го порядка применительно к решению уравнения Боука — Вена
можно записать в векторном виде:
Γj+1 = Γj +
1
6
(k1
j + 2k2
j + 2k3
j + k4
j), (61)
где
k1
j = htFj|
ξj
i=ξj
i
,
k2
j = htFj|
ξj
i=ξj
i+k1
j
i/2,
k3
j = htFj|
ξj
i=ξj
i+k2
j
i/2,
k4
j = htFj|
ξj
i=ξj
i+k3
j
i
,
где Fj, k1
j, k2
j, k3
j, k4
j есть векторы-столбцы на j временном интервале. Таким образом,
представленный метод позволяет оптимизировать вычисления, значительно сократив время их
выполнения (в средах, где есть поддержка векторного исчисления, например, среда Mathlab
или библиотека NumPy для языка Python и др.).
4. Результаты моделирования колебаний
В настоящем разделе описываются результаты моделирования колебаний классической
и гистерезисной балок в зависимости от разных типов нагрузок. Моделирование феномена
гистерезиса осуществляется посредством двух преобразователей: Боука — Вена и Прандтля —
Ишлинского. Также в разделе приводится сравнительная характеристика между этими моде-
лями. Начальные условия полагаются нулевыми (балка покоится в начальный момент времени
и не имеет начальной скорости соответственно) и ненулевыми граничные условия, когда моде-
лируется сейсмическая нагрузка (на концах моделируется воздействие сейсмической волны) и
нулевыми для иных нагрузок (рис. 9). Характеристики балки постоянны во всех экспериментах
(кроме тех случаев, где явно указаны другие параметры):
• длина балки балки равна L = 1;
• время моделирования равно T = 5;
• коэффициент a2 равняется 0.5;
параметры численного моделирования полагались следующими:
• шаг по времени равен ht = 0.00125;
Карпов Е.А. Динамика балки Эйлера — Бернулли с учетом распределенных гистерезисных свойств
Karpov E.A. Dynamics of the Euler — Bernoully beam with distributed hysteresis properties 50 из 62
• шаг по пространсву равен 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)
Seismic S(t) influence: S(t) S(t)
Рис. 9. Шарнирное закрепление балки с учетом внешней нагрузки 1(x, t) (управления)
и сейсмологической нагрузки S(t) в месте закрепления
Fig. 9. Fixed boundary of the beam taking into account the external load 1(x, t) (control)
and seismological load S(t) at the boundary connection
4.1. Моделирование нагрузки как функции периодической структуры
Одним из классических типов нагрузки является периодическая нагрузка. С ее помощью
можно понять, как реагирует система на периодические воздействия с меняющейся частотой
или амплитудой. Обычно в качестве периодической нагрузки принято использовать триго-
нометрические функции, изменяющиеся в зависимости от времени, вида A sin(ωt + ϕ0) или
A cos(ωt + ϕ0), где A — амплитуда, ω — частота и ϕ0 — начальная фаза колебаний. В настоя-
щей статье в качестве функции периодического воздействия будет использоваться функция
1(x, t) = 0.12 sin(2π/2 − 4πt) (рис. 10).
Как можно видеть из полученных результатов (рис. 11), реакция каждого типа балки
на периодическое воздействие разная. Классическая балка с начального момента времени
начинает демонстрировать квазипериодические формы колебаний c двумя частотами (основной
с наибольшей мощностью, и второстепенная, с меньшей), что подтверждается анализом спектра
Фурье-преобразования (рис. 12, а). Балка с гистерезисными свойствами, формализуемыми
преобразователем Боука — Вена, демонстрирует схожие результаты с классической балкой.
А именно на начальном временном отрезке можно видеть схожие колебания с классической
моделью. Однако, начиная с некоторого момента времени, устанавливаются периодические коле-
бания с постоянной частотой (отметим, что частота у каждого из преобразователей отличается,
ровно как и амплитуда колебаний). Анализ спектра Фурье-преобразования показывает, что
использование преобразователя Боука — Вена способствует уменьшению мощности неосновной
частоты колебаний (рис. 12, в) практически нивелируя ее. Используя в качестве носителя
гистерезисных свойств преобразователь Прандтля — Ишлинского, можно добиться лучших
Вестник Самарского университета. Естественнонаучная серия 2024. Том 3, № 30. С. 35–62
Vestnik of Samara State University. Natural Science Series 2024, vol. 3, no. 30, pp. 35–62 51 из 62
Рис. 10. График функции нагрузки 1(x, t) = 0.12sin(π/2 − 4πt)
Fig. 10. Graphical representation of the load function 1(x, t) = 0.12sin(π/2 − 4πt)

 
 

-



-



-






















(/ )
а

 
 

-



-



-






















(/ )
б
-


-


-
















-


-



-





















(/     
)



(
/ )
в
-

 -









-


-












[



] / )



(
/ )
г
Рис. 11. Графики колебаний балки в точке x = L/2 в зависимости от времени: синяя кривая —
классическая модель балки, оранжевая кривая — гистерезисная балка с использованием
преобразователя Прандтля — Ишлинского, красная кривая — гистерезисная балка с использованием
преобразователя Боука — Вена (а, б); параметрические графики гистерезисных петель (время t как
параметр): входно-выходная зависимость входа uxx от выхода Γ[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 Γ[uxx] at the point x = L/2 (c, d)
результатов в сравнении с остальными моделями. А именно с ее помощью, во-первых, можно
добиться уменьшения амплитуды основной частоты, во-вторых, редуцировать второстепенную
частоту (рис. 12, б). Представленные гистерезисные зависимости (рис. 11, в, г) позволяют
Карпов Е.А. Динамика балки Эйлера — Бернулли с учетом распределенных гистерезисных свойств
Karpov E.A. Dynamics of the Euler — Bernoully beam with distributed hysteresis properties 52 из 62
оценить энергетический вклад гистерезисного звена в редуцирование колебаний: площадь петли
пропорциональна энергии, “отводимой” гистерезисным звеном из системы.
Таким образом, преобразователь Прандтля — Ишлинского позволяет избавить колебания
балки от квазипериодической структуры, трансформируя их в чисто периодическую, а преоб-
разователь Боука — Вена способствует уменьшению мощности неосновной частоты колебаний
балки.


































    
  
   
    
а










































    
       

б































   




       
в
Рис. 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. Моделирование воздействия импульсной нагрузки
Нередко в рамках прикладных задач важно знать реакцию несущей конструкции на им-
пульсную нагрузку. Самым простым примером такой нагрузки является точечный удар по
конструкции с некоторой “силой”. Классическим методом формализации импульсной нагрузки
является дельта-функция Дирака. Однако, принимая во внимание использование численных
методов в работе, имеет смысл рассмотреть подобную ей функцию, которая в одномерном
случае имеет вид:
ˆ δ(t) =

A, если t = 0,
0, иначе,
а в двумерном:
ˆ δ(x, t) =

A, если t = 0 и x = 0,
0, иначе,
где A — “сила” удара. Использование такой функции позволяет нивелировать бесконечности в
расчетах, что сильно упрощает задачу моделирования. Ясно, что балка за счет внутренних сил
должна “поглощать” удар, преобразуя энергию удара в иные ее формы, например, в тепловую
энергею (нагрев балки). В настоящей статье исследуется реакция классической и “гистерезисной”
балок на импульсную нагрузку (удар), приложенную к ее определенной дискретной точке в
некоторый момент времени. Таким образом, функцию 1(x, t) можно представить в виде:
1(x, t) = ˆ δ(x − χ, t − τ) = ˆ δ(x − χ) ˆ δ(t − τ),
где χ ∈ [0, L] есть дискретная точка балки, в которой действует импульсная нагрузка, а
τ ∈ (0, T) — момент времени, в который действует нагрузка.
Результаты моделирования, продемонстрированные на рис. 13, показали, что учет гисте-
резисных свойств приводит к совершенно иной динамике колебаний, нежели в классической
модели. А именно начиная с момента времени t = τ классическая балка начинает совершать
колебания с постоянной частотой, причем в ходе экспериментов было установлено, что им-
пульсная нагрузка любой “силы” выводит балку из состояния равновесия, заставляя балку
Вестник Самарского университета. Естественнонаучная серия 2024. Том 3, № 30. С. 35–62
Vestnik of Samara State University. Natural Science Series 2024, vol. 3, no. 30, pp. 35–62 53 из 62

 
 

-


-




















(/ )
а

 
 

-


-























(/ )
б

 
 



















(/ )
в
Рис. 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
бесконечно долго колебаться. Напротив, балка с учетом гистерезисных свойств с момента
времени t = τ начинает отклоняться от положения равновесия, достигает некоторой точки,
совершая при этом одиночное колебание малой амплитуды, после чего достигает устойчивого
положения, отличного от положения равновесия и остается в этом положении бесконечно долго.
При этом на рис. 13, в можно видеть, что результаты, демонстрируемые с использованием
преобразователей Прандтля — Ишлинского и Боука — Вена, отличаются, однако схожи по
своей структуре. Схожесть результатов подтверждается и подобностью гистерезисных петель,
продемонстрированных на рис. 14.
-

-





-
-

-






(/     )



(
/ )
а
-
 -

-

 -

 -

 -





-

-

-

-

-






(    /
  )



(
/ )
б
Рис. 14. Параметрические графики гистерезисных петель (время t как параметр): входно-выходная
зависимость входа uxx от выхода Γ[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 Γ[uxx] at the point x = L/2: a — the Prandtl — Ishlinskii operator, b — the Bouc — Wen
operator
Наблюдаемое поведение колебаний у балки с учетом гистерезисных свойств можно интер-
претировать как более удачное, чем у балки, моделируемой классическим методом. Такой
вывод можно сделать, полагая, что, получая удар, балка деформируется в сторону удара
(соседние точки деформируются, в том числе образуя адекватный прогиб, см. рис. 15), при
этом гистерезис способствует диссипации энергии удара, нивелируя дальнейшие колебания
балки, в отличие от классической модели, колебания которой имеют периодический характер.
Карпов Е.А. Динамика балки Эйлера — Бернулли с учетом распределенных гистерезисных свойств
Karpov E.A. Dynamics of the Euler — Bernoully beam with distributed hysteresis properties 54 из 62






 




























( /)
Рис. 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”, имеющая вид:
Ψa,b(t) = Ψ
 
t − b
a
!
, (62)
где a и b являются параметрами растяжения и сдвига материнской функции, соответственно.
Под вейвлет-функциями в рамках настоящей работы понимается хорошо известное семейство
вейвлетов Добеши.
Начнем анализ полученных результатов со сравнения классической и “гистерезисной” балок
с гистерезисным преобразователем Боука — Вена. Как видно из рис. 16 на начальных вре-
менных интервалах (с 0 по 3 отчет модельного времени) балки ведут себя одинаково, однако
спустя 8 интервала временного отчета “гистерезисная” балка стаблизируется, в то время как
классическая балка демонстрирует периодические колебания. Такая структура колебаний
объясняется диссипирующими свойствами гистерезисного преобразователя (рис. 17), который
трансформирует энергию колебаний сейсмической волны в иные формы.
Также из результатов моделирования, представленных на рис. 16, видно, что реация клас-
сической и “гистерезисной” балок различна. А именно в моментах воздействия сейсмической
волны гистерезисная балка демонстрирует возмущения меньшой амплитуды (особенно четко это
наблюдается для балки с использованием модели Прандтля — Ишлинского) в сравнении с клас-
сической балкой. После прекращения сейсмического воздействия обе балки стабилизируются в
состояние равновесия, однако гистерезисная балка, формализуемая преобразователем Прандт-
ля — Ишлинского, оказывается в положении равновесия, отличном от исходного (рис. 16, б),
что согласуется с упруго-пластическими свойствами гистерезисного материала.
Также, как видно из представленных результатов, обе балки все еще продолжают совер-
шать колебания после воздействия нагрузки. Однако балка с учетом гистерезисных свойств
показывает колебания значительно меньшей амплитуды, чем классическая (0.0004 против 0.1
соответственно).
Из продемонстрированных результатов можно сделать вывод, что учет гистерезисных
свойств позволяет идентифицировать новые формы поведения в сравнении с классической
моделью.
Вестник Самарского университета. Естественнонаучная серия 2024. Том 3, № 30. С. 35–62
Vestnik of Samara State University. Natural Science Series 2024, vol. 3, no. 30, pp. 35–62 55 из 62

 
 

-

-









а

 
 

-

-










(/ )
б
Рис. 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 как параметр): входно-выходная зависимость входа uxx от выхода Γ[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 Γ[uxx] at the point x = L/2: a — the Prandtl — Ishlinskii operator, b — the Bouc — Wen operator
Заключение
В настоящей статье исследована динамика балки с распределенными гистерезисными свойствами. В первой части работы получены уравнения с использованием гамильтоновского подхода, описывающие модель движений балки, когда носитель гистерезисных свойств распределен по всей ее длине. Приведены результаты моделирования колебаний балки, находящейся под внешней нагрузкой различной природы: сейсмической (формализуемой материнской вейвлет-функцией из семейства вейвлетов Добеши), ударной и периодической. Произведено сравнение характера колебаний классической и “гистерезисной” балок. Для описания гистерезисных свойств использовались два подхода: феноменологический (преобразователь Боука — Вена) и конструктивный (преобразователь Прандтля — Ишлинского). На основе численных экспериментов установлено, что балка с распределенными гистерезисными свойствами обладает повышенной (по сравнению с классической моделью) устойчивостью по отношению к внешним воздействиям различной природы. Иными словами, области устойчивости (в пространстве
параметров этой системы) шире для гистерезисной балки. Полученные результаты могут найти применение при разработке и проектировании несущих конструкций в сейсмоопасных регионах