An application of Mueller's method for determining eigenfrequencies of vibrations of viscoelastic bodies with frequency-dependent characteristics of a material

Full Text

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

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

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

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

Демпфирующие свойства конструкции можно определить, основываясь на учете временного фактора в рамках теорий сплошной среды наследственного вида: вязкого сопротивления Максвелла, вязкого трения Кельвина – Фойгта, наследственности Больцмана – Вольтерра, реологической модели Пойтинга – Томсона, термодиффузионной теории Зинера и других [2]. Формулировка общего вида и конкретизация соотношений для изотропных, анизотропных и неоднородных материалов, развитие общих методов вязкоупругости содержатся в фундаментальных монографиях А. А. Ильюшина и Б. Е. Победри [2], Ю. Н. Работнова [3, 4], Н. Х. Арутюняна [5, 6], И. И. Бугакова [7], В. Г. Карнаухова [8], М. А. Колтунова [9], В. В. Москвитина [10], А. Р. Ржаницына [11, 12], Д. Бленда [13], Р. Кристенсена [14], Дж. Ферри [15] и других. Моделям, содержащим дробные производные и другие операторы дробного порядка, посвящены обзоры [16] и [17], в которых рассмотрено 337 и 250 работ соответственно.

Наиболее типичными представителями такого рода материалов являются полимерные и композиционные материалы, горные породы, бетонные и железобетонные конструкции. Здесь использовалась наиболее общая линейная теория Больцмана – Вольтерра, отражающая практически все особенности квазистатического и динамического поведения вязкоупругих материалов [18–20]. Использование соотношений данной теории в задачах о собственных или вынужденных установившихся колебаниях позволяет осуществить переход к описанию материальных характеристик в терминах комплексных динамических модулей. В этом случае определяющие соотношения будут иметь точно такой же вид, как и для упругого материала, но при этом коэффициенты, определяющие связь между компонентами тензоров напряжений и деформаций, будут являться комплексными величинами.

Решением задачи о собственных колебаниях кусочно-однородных вязкоупругих тел в приведенной постановке являются комплексные собственные частоты колебаний, действительная часть которых определяет частоту, а мнимая — показатель демпфирования (скорость затухания).

Использование комплексных динамических модулей в задачах о вынужденных установившихся колебаниях существенно упрощает процесс получения их решений как аналитическими [21–23], так и численными методами, в частности методом конечных элементов [21, 22, 24, 25]. Более того, используемые в данном случае комплексные величины модулей могут быть получены из экспериментов напрямую [26], что является большим преимуществом данного подхода при использовании в практических приложениях. На сегодняшний день процедуры численного решения задач о вынужденных установившихся колебаниях реализованы во многих пакетах конечно-элементного анализа, в том числе в ANSYS.

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

Как отмечается в работе [27], на примере материала типа «ПДИ», который является механическим аналогом (имитатором) твердого ракетного топлива и представляет собой высоконаполненную полимерную композицию [10], экспериментально показано, что динамический модуль является монотонно возрастающей функцией логарифма частоты нагружения и в частотном диапазоне от 0 до 200 Гц с увеличением частоты увеличивается более чем в 5 раз.

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

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

Данным условиям наилучшим образом удовлетворяет метод Мюллера, который является основой представленного в данной работе способа решения задачи о собственных колебаниях вязкоупругих конструкций с учетом частотной зависимости комплексных динамических модулей, описывающих вязкоупругие свойства материала конструкции. Метод Мюллера, варианты выбора начальных приближений для решения алгебраической проблемы действительных или комплексных собственных значений подробно рассмотрены в работах [33, 34].

Метод Мюллера зарекомендовал себя как весьма универсальный и подходящий для решения широкого класса спектральных задач, сводящихся к алгебраическим проблемам на комплексные собственные значения [35–40]. Более того, стоит отметить, что данный метод позволяет получить решения задач, в характеристические уравнения которых искомый параметр входит сложным образом (различного рода степенные зависимости, обратная зависимость и т.д.) [38–40]. Однако несмотря на все описанные выше особенности, использовать метод Мюллера напрямую для решения задачи о собственных колебаниях с учетом зависимости механических характеристик вязкоупругого материала от частоты оказалось невозможно.

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

В данной работе представлен вариант алгоритма определения комплексных собственных частот колебаний кусочно-однородных вязкоупругих тел, выполненных из материалов с частотно-зависимыми характеристиками. Он является дальнейшим развитием подходов к построению численных решений задачи о собственных колебаниях кусочно-однородных упругих или вязкоупругих тел, развиваемых под руководством В. П. Матвеенко [19, 20]. Работоспособность и эффективность предложенного метода продемонстрированы на примере двухслойной консольно защемленной пластинки, один слой которой выполнен из упругого материала, а другой — из вязкоупругого. Достоверность полученных результатов подтверждается сравнением собственных частот колебаний, определенных решением задачи о собственных колебаниях такого рода конструкций, с резонансными частотами на амплитудно-частотных характеристиках перемещений, полученных решением задачи об установившихся вынужденных колебаниях в пакете прикладных программ ANSYS.

1. Математическая постановка задач о собственных и вынужденных установившихся колебаниях кусочно-однородных вязкоупругих тел. Особенности математической постановки задачи о собственных колебаниях кусочно-однородных тел, имеющих в своем составе элементы, выполненные из вязкоупругого материала, механическое поведение которого описывается линейной наследственной теорией вязкоупругости с учетом зависимости механических характеристик от частоты внешнего воздействия, приводят к необходимости разработки специальных подходов к построению численных решений на основе метода конечных элементов. Для этого воспользуемся способом, сформулированным ранее в работах [18–20].

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

$\sum _{n=1}^N{\int }_{V_1^n}({\sigma }_{ij}\delta {\epsilon }_{ij}+{\rho }_n{\ddot{u}}_i\delta u_i-f_i\delta u_i)dV+$$\sum _{m=1}^M{\int }_{V_2^m}({\sigma }_{ij}\delta {\epsilon }_{ij}+{\rho }_m{\ddot{u}}_i\delta u_i-f_i\delta u_i)dV-{\int }_{S_{\sigma }}{p}_i\delta u_{i}dS=0$ (1)

В случае кусочно-однородного тела объемом $V{V}_1\cup V_2$ часть объема $V_1$ состоит из$N$ упругих элементов, а часть объема $V_2$ — из $M$вязкоупругих элементов. В уравнение (1) входят компоненты: ${\sigma }_{ij}$ — симметричного тензора напряжений Коши; ${\epsilon }_{ij}\frac{1}{2}(u_{i\mathrm{,}j}+u_{j\mathrm{,}i})$ — симметричного тензора деформаций; $u_i$ — вектора перемещений; $f_i$ — вектора объемных усилий; вектора поверхностных усилий $p_i$, заданных на $S_{\sigma }$ — части поверхности кусочно-однородного тела объемом $V$; ${\rho }_n$, ${\rho }_m$ — удельные плотности упругих и вязкоупругих частей кусочно-однородного тела.

Упругие и вязкоупругие слои идеально скреплены между собой. На поверхности, ограничивающей объем рассматриваемого кусочно-однородного тела, заданы кинематические и силовые граничные условия:

$u=U_i на S_u\mathrm{,} {\sigma }_{ij}{n}_j=p_i на S_{\sigma }\mathrm{,}$ (2)

где $n_j$ — компоненты вектора внешней нормали к поверхности $S_{\sigma }$, $S_u$ — часть поверхности кусочно-однородного тела объемом $V$, на которой заданы кинематические граничные условия.

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

    • при рассмотрении собственных колебаний компоненты вектора поверхностных усилий и вектора перемещений принимают нулевые значения: $U_i$= 0, $p_i = 0$;

• в случае вынужденных установившихся колебаний компоненты векторов перемещений и поверхностных усилий изменяются по гармоническому закону:

$U_i=U_i^0\cos \Omega t$, $p_i=p_i^0\cos \Omega t$.

Для упругих составляющих объема $V_1$ при малых деформациях выполняются физические соотношения для изотропного [2, 41]

${\sigma }_{ij}-\sigma {\delta }_{ij}=2G^{\left(n\right)}({\epsilon }_{ij}-\frac{1}{3}\theta {\delta }_{ij}), \sigma =B^{\left(n\right)}\theta$ (3) 

или анизотропного материалов [2, 41]

${\sigma }_{ij}=C_{ijkl}^{\left(n\right)}{\epsilon }_{kl}\mathrm{.}$ (4)

Здесь $G^{\left(n\right)}$, $B^{\left(n\right)}$ — упругие сдвиговые и объемные модули, $\sigma$ — среднее напряжение, $\theta$ — объемная деформация, $C_{ijkl}^{\left(n\right)}$ — компоненты тензора упругих констант анизотропного тела.

Для вязкоупругих частей рассматриваемого тела выполняются соотношения линейной вязкоупругости для изотропного материала [2–20]:

$\begin{array}{l}{\sigma }_{ij}-\sigma {\delta }_{ij}=2G_0^{\left(m\right)}(e_{ij}-{\int }_{-\infty }^t{H}^{\left(m\right)}(t-\tau \left)e_{ij}\right(\tau \left)d\tau \right),\\ \sigma B_0^{\left(m\right)}(\theta -{\int }_{-\infty }^t{F}^{\left(m\right)}(t-\tau \left)\theta \right(\tau \left)d\tau \right).\end{array}$(5)

Здесь $G_0^{\left(m\right)}$, $B_0^{\left(m\right)}$ — мгновенные сдвиговые и объемные модули; $s_{ij}$, $e_{ij}$ — компоненты девиаторов тензоров напряжений и деформаций; $F^{\left(m\right)}$, $H^{\left(m\right)}$ — ядра релаксации.

В рамках данной работы будут рассматриваться как собственные, так и вынужденные колебания кусочно-однородного вязкоупругого тела. В случае вынужденных колебаний решение задачи, описываемой уравнением (1), будет отыскиваться в виде

$\overline{u}\left(\overline{x}\mathrm{,}t\right)={\overline{u}}_0\left(\overline{x}\right)e^{i\Omega t}$. (6)

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

$\overline{u}\left(\overline{x}\mathrm{,}t\right)={\overline{u}}_0\left(\overline{x}\right)e^{i\omega t}\mathrm{.}$ (7)

В выражениях (6), (7) $\overline{u}\left(\overline{x}\mathrm{,}t\right)$ — обобщенный вектор состояния; ${\overline{u}}_0\left(\overline{x}\right)$ — вектор перемещений, зависящий только от пространственных координат (форма колебаний); $\Omega$ — частота внешнего возбуждения; $\omega ={\omega }_{Re}+i{\omega }_{Im}$ комплексная собственная частота колебаний, в которой действительная часть ${\omega }_{Re}$ является круговой собственной частотой колебаний, а мнимая часть ${\omega }_{Im}$ характеризует скорость их затухания.

Переход к описанию вязкоупругих свойств материала в терминах комплексных динамических модулей в задачах о собственных и вынужденных установившихся колебаниях требует некоторых математических преобразований. При вынужденных установившихся колебаниях нижний предел интегрирования в соотношениях (5) заменяется на $-\infty$ [20]. После подстановки решения (7) в (5), замены переменной $s=t-\tau$ ($\tau =t-s,$ $d\tau =ds$) и соответствующих математических преобразований уравнения, определяющие связь между компонентами тензоров напряжений и деформаций, в случае изотропного материала принимают следующий вид:

$\begin{array}{l}{s}_{ij}^0=2G_0^{\left(m\right)}{e}_{ij}^0\left(x_i\right)(1-{\int }_{-\infty }^t{H}^{\left(m\right)}(s)\cos \Omega s ds+i{\int }_{-\infty }^t{H}^{\left(m\right)}(s)\sin \Omega s ds)\\ {\sigma }^0=B_0^{\left(m\right)}{\theta }^0\left(x_i\right)(1-{\int }_{-\infty }^t{F}^{\left(m\right)}(s)\cos \Omega s ds+i{\int }_{-\infty }^t{F}^{\left(m\right)}(s)\sin \Omega s ds).\end{array}$

Далее вводятся следующие обозначения, которые определяют действительную и мнимую компоненты комплексных динамических модулей сдвига и объемного сжатия [18–20]:

$G_{Re}^{\left(m\right)}=G_0^{\left(m\right)}(1-{\int }_{-\infty }^t{H}^{\left(m\right)}(s)\cos \Omega s ds),$

$G_{Im}^{\left(m\right)}=G_0^{\left(m\right)}{\int }_{-\infty }^t{H}^{\left(m\right)}\left(s\right)\sin \Omega s ds\mathrm{,}$ (8)

$B_{Re}^{\left(m\right)}=B_0^{\left(m\right)}\left(1-{\int }_{-\infty }^t{F}^{\left(m\right)}\left(s\right)\cos \Omega s ds\right)\mathrm{,}$

$B_{Im}^{\left(m\right)}=B_0^{\left(m\right)}{\int }_{-\infty }^t{F}^{\left(m\right)}\left(s\right)\sin \Omega s ds\mathrm{.}$

С учетом (8) определяющие соотношения (5) для вынужденных установившихся колебаний можно переписать в следующем виде:

${\overline{s}}_{ij}=2{\tilde{G}}^{\left(m\right)}{e}_{ij}\mathrm{,} \overline{\sigma }={\tilde{B}}^{\left(m\right)}\theta \mathrm{,}$ (9)

где ${\tilde{G}}^{\left(m\right)}=G_{Re}^{\left(m\right)}+i{G}_{Im}^{\left(m\right)},$ ${\tilde{B}}^{\left(m\right)}=B_{Re}^{\left(m\right)}+i{B}_{Im}^{\left(m\right)}$.

Аналогичные математические преобразования можно осуществить и для задачи о собственных колебаниях. Для этого необходимо подставить в (7) выражение для комплексной собственной частоты $\omega ={\omega }_{Re}+i{\omega }_{Im}$. В результате с учетом этой подстановки искомое решение задачи о собственных колебаниях примет следующий вид:

$\overline{u}\left(\overline{x}\mathrm{,}t\right)={\overline{u}}_0\left(\overline{x}\right)e^{i{\omega }_{Re}t}{e}^{-{\omega }_{Im}t}\mathrm{.}$

Однако для получения конечного вида определяющих соотношений в случае собственных колебаний вводятся два допущения [18–20]:

• колебания происходят с медленно меняющимися амплитудами;
• начальные возмущения не влияют на поведение системы.

Первое допущение позволяет вынести из-под знака интеграла $e^{-{\omega }_{Im}t}$. Второе допущение позволяет изменить нижний предел интегрирования на $-\infty$. Подробное обоснование сделанных предположений содержится в работах [19, 20].

Осуществив математические преобразования, аналогичные тем, что были проделаны для случая вынужденных установившихся колебаний, мы получаем физические уравнения в виде 

$s_{ij}^0=2G_0^{\left(m\right)}{e}_{ij}^0\left(x_i\right)(1-{\int }_{-\infty }^t{H}^{\left(m\right)}(s)\cos {\omega }_{Re}s ds+i{\int }_{-\infty }^t{H}^{\left(m\right)}(s)\sin {\omega }_{Re}s ds),$

${\sigma }^0=B_0^{\left(m\right)}{\theta }^0\left(x_i\right)(1-{\int }_{-\infty }^t{F}^{m}s\cos {\omega }_{Re}s ds+i{\int }_{-\infty }^t{F}^{\left(m\right)}(s)\sin {\omega }_{Re}s ds).$

Далее вводятся обозначения для компонент комплексных модулей, аналогичные (8):

$G_{Re}^{\left(m\right)}=G_0^{\left(m\right)}(1-{\int }_{-\infty }^t{H}^{\left(m\right)}(s)\cos {\omega }_{Re}s ds),$

$G_{Im}^{\left(m\right)}=G_0^{\left(m\right)}{\int }_{-\infty }^t{H}^{\left(m\right)}\left(s\right)\sin {\omega }_{Re}s ds\mathrm{,}$(10)

$B_{Re}^{\left(m\right)}=B_0^{\left(m\right)}(1-{\int }_{-\infty }^t{F}^{\left(m\right)}(s)\cos {\omega }_{Re}s ds),$

$B_{Im}^{\left(m\right)}=B_0^{\left(m\right)}{\int }_{-\infty }^t{F}^{\left(m\right)}\left(s\right)\sin {\omega }_{Re}s ds\mathrm{,}$

которые определяют действительную и мнимую части комплексного динамического модуля. Соотношения (4) для задачи о собственных колебаниях примут вид

${\overline{s}}_{ij}\approx 2{\tilde{G}}^{\left(m\right)}{e}_{ij}\mathrm{,} \overline{\sigma }\approx {\tilde{B}}^{\left(m\right)}\theta \mathrm{,} \left(11\right)$

где ${\tilde{G}}^{\left(m\right)}=G_{Re}^{\left(m\right)}+i{G}_{Im}^{\left(m\right)},$${\tilde{B}}^{\left(m\right)}=B_{Re}^{\left(m\right)}+i{B}_{Im}^{\left(m\right)}.$

Здесь следует отметить, что соотношения (9) выполняются точно, а (11) — приближенно. Анализ соотношений (8) и (10) позволяет сделать вывод о том, что с точки зрения физического смысла данные соотношения эквивалентны. Тем не менее данные соотношения имеют существенное отличие, состоящее в том, что в общем случае выражения (8) являются функциями вещественного параметра $\Omega$, а выражения (10) — функциями действительной части ${\omega }_{Re}$комплексной собственной частоты $\omega ={\omega }_{Re}+i{\omega }_{Im}$. Данная особенность приводит к тому, что процесс построения алгоритма численного решения задач о собственных колебаниях методом конечных элементов (МКЭ) оказывается крайне нетривиальной задачей.

2. Численная реализация задач о собственных и о вынужденных установившихся колебаниях вязкоупругих тел с частотно-зависимыми модулями методом конечных элементов. Решение сформулированных выше вариационных задач о собственных и о вынужденных установившихся колебаниях вязкоупругих тел с частотно-зависимыми комплексными динамическими модулями может быть осуществлено методом конечных элементов (МКЭ). Типовые процедуры МКЭ приводят рассматриваемые задачи к системам линейных алгебраических уравнений, которые могут быть записаны в матричной форме (12) или (13) для соответствующих задач: 

$\left(\left[K\right]-{\omega }^2\left[M\right]\right)\left\{\xi \right\}=\left\{0\right\}, \left(12\right)$

$\left(\left[K\right]-{\Omega }^2\left[M\right]\right)\left\{\xi \right\}=\left\{F\right\}, \left(13\right)$

где $\left[K\right]$ — глобальная матрица жесткости, $\left[M\right]$ — глобальная матрица масс, $\left\{\xi \right\}$ — вектор узловых неизвестных, $\left\{F\right\}$ — глобальный вектор узловых нагрузок.

В случае кусочно-однородного вязкоупругого тела, состоящего из упругих и вязкоупругих частей, глобальную матрицу жесткости в уравнениях (12)–(13) можно представить в виде суммы упругой $\left[K_{elast}\right]$ и вязкоупругой $\left[K_{vis}\right]$ составляющих:

$\left[K\right]=\left[K_{vis}\right]+\left[K_{elast}\right].$

В свою очередь, в случае изотропного материала для вязкоупругих компонент рассматриваемого кусочно-однородного тела матрицу жесткости вязкоупругой части, согласно соотношениям (8), (9) и (10), (11), можно представить в виде (14), воспользовавшись алгоритмом, приведенным в [35]:

$\left[K_{vis}\right]=\widetilde{B[}{K}_B]+\widetilde{G[}{K}_G], \left(14\right)$

в котором комплексные динамические модули $\tilde{G}$, $\tilde{B}$ в общем случае являются либо функциями частоты внешнего возбуждения $\Omega$ (при рассмотрении вынужденных установившихся колебаний), либо функциями действительной части комплексной собственной частоты ${\omega }_{Re}$ (в случае анализа собственных колебаний). С учетом введенных обозначений матричные уравнения (12)–(13) можно переписать в следующем виде:

$\left(\right[K_{elast}]+\widetilde{B[}{K}_B]\left({\omega }_{Re}\right)+\widetilde{G[}{K}_G\left]\right({\omega }_{Re})-{\omega }^2[M\left]\right)\left\{\xi \right\}=\left\{0\right\}, \left(15\right)$

$\left(\right[K_{elast}]+\tilde{B}[K_B\left]\right(\Omega )+\tilde{G}[K_G\left]\right(\Omega )-{\Omega }^2[M\left]\right)\left\{\xi \right\}=\left\{F\right\}. \left(16\right)$

Здесь выражение $\widetilde{G[}{K}_G]+\widetilde{B[}{K}_B]$ описывает вязкоупругую часть рассматриваемого кусочно-однородного тела объемом $V_2$, $\left[K_B\right]$, $\left[K_G\right]$ — объемная и сдвиговая компоненты матрицы жесткости вязкоупругой части.

При рассмотрении задачи о вынужденных установившихся колебаниях, описываемой матричным уравнением (16), ее решение в виде значений составляющих вектора узловых неизвестных $\left\{\xi \right\}$ для любого заданного значения частоты колебаний $\Omega$ отыскивается путем решения системы линейных алгебраических уравнений (СЛАУ) при помощи любого известного метода. В частности, для решения данной задачи могут быть использованы коммерческие пакеты прикладных программ, например ANSYS.

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

Однако в случае рассмотрения собственных колебаний получить решение задачи, описываемой матричным уравнением (15), в виде значений комплексных собственных частот $\omega ={\omega }_{Re}+i{\omega }_{Im}$ и форм колебаний $\left\{\xi \right\}=\left\{{\xi }_{Re}\right\}+i\left\{{\xi }_{Im}\right\}$ оказывается весьма затруднительно. Согласно определению, условием существования нетривиального решения $\left\{\xi \right\}\ne \left\{0\right\}$ для СЛАУ, описываемой матричным уравнением (15), является равенство нулю определителя матрицы коэффициентов системы:

$D=\left(\right[K_{elast}]+\tilde{B}[K_B]+\tilde{G}[K_G]-{\omega }^2[M\left]\right)=0.$

Таким образом, мы получили алгебраическую проблему комплексных собственных значений для матрицы $\widetilde{\left[G\right[}{K}_G]+\tilde{B}[K_B]+[K_{elast}]-{\omega }^2[M\left]\right].$

При выборе метода решения алгебраической проблемы собственных значений необходимо учитывать большую размерность СЛАУ и возможность решения частичной алгебраической проблемы комплексных собственных значений. Следует отметить, что для рассматриваемой задачи, описываемой матричным уравнением (15), имеет смысл решать только частичную проблему собственных значений, то есть определять не весь спектр собственных частот, а только какую-либо его часть. Последнее обстоятельство выдвигает жесткое условие, состоящее в том, что собственные значения должны быть гарантированно определены либо в заданном интервале, либо в требуемой последовательности, например, в порядке возрастания действительных частей в случае комплексных собственных значений.

Исследования динамического поведения конструкций в настоящее время активно проводятся с помощью алгоритмов, в той или иной мере реализованных в коммерческих пакетах программ инженерных расчетов, реализующих метод конечных элементов (MatLab, ANSYS, Abaqus, MSC NASTRAN и т.д.). Однако решение задачи о собственных колебаниях вязкоупругих конструкций, материал которых описывается наследственной теорией вязкоупругости на основе комплексных динамических модулей, изменяющихся независимо друг от друга, стандартными методами в пакете ANSYS невозможно.

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

Для случая, когда составляющие комплексных динамических модулей сдвига и объемного сжатия не зависят от частоты колебаний, т.е. являются постоянными $\tilde{G}=const$, $\tilde{B}=const$ (например, на ограниченном участке рассматриваемого частотного диапазона), разработан алгоритм численного решения задачи о собственных колебаниях, построенный на использовании возможностей пакета ANSYS, которые дают возможность записывать глобальные ансамблированные матрицы жесткости и масс, позволяющие сформировать уравнение вида (15).

Для решения алгебраической проблемы собственных значений для матриц с постоянными коэффициентами в рамках данного алгоритма используется метод Мюллера с различными сценариями выбора начальных приближений [33, 34]. В результате решения проблемы собственных значений методом Мюллера получается спектр комплексных собственных частот колебаний вязкоупругой конструкции с частотно-независимыми характеристиками материала. Так как в общем случае характеристики вязкоупругого материала зависят от частоты колебаний внешнего воздействия, коэффициенты матрицы жесткости также являются переменными величинами функциями частоты, и напрямую воспользоваться каким-либо методом для решения частичной алгебраической проблемы комплексных собственных значений (16) не представляется возможным. Авторам не удалось найти в литературе сведений о способах решения данной проблемы для наиболее общих и широко применяемых для большого класса материалов моделей вязкоупругости. При этом необходимо отметить, что предлагаемый авторами алгоритм является универсальным и подходящим для решения модальных задач при использовании различных моделей вязкоупругости, в том числе и моделей, содержащих дробные производные и другие операторы дробного порядка [16, 17].

В научной литературе [21, 22, 26–32, 43, 44] приводятся различные варианты графических зависимостей составляющих комплексного динамического модуля от частоты колебаний. Традиционно такого рода зависимости модулей сдвига и объемного сжатия определяются экспериментально на основе метода динамического механического анализа (ДМА), и строятся они в зависимости от частоты внешнего возбуждения. При переходе к задаче о собственных колебаниях данные зависимости будут являться уже функциями действительной части комплексной собственной частоты ${\omega }_{Re}$ (15), которая имеет смысл круговой собственной частоты колебаний. Однако ${\omega }_{Re}$ является лишь компонентой искомой комплексной собственной частоты $\omega ={\omega }_{Re}+i{\omega }_{Im}$, что приводит к дополнительным сложностям при реализации численного решения рассматриваемой задачи. Очевидно, что в общем случае собственные значения матрицы $\left[\widetilde{G(}{\omega }_{Re}\right)\left[K_G\right]+\tilde{B}\left({\omega }_{Re}\right)\left[K_B\right]+\left[K_{elast}\right]-{\omega }^2\left[M\right]]$

 будут отличаться от собственных значений матрицы

$\left[\widetilde{G(}\omega \right)\left[K_G\right]+\widetilde{B(}\omega \left)\right[K_B]+[K_{elast}]-{\omega }^2[M\left]\right].$

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

3. Алгоритм определения значений собственных частот колебаний кусочно-однородных вязкоупругих тел, выполненных из материалов с частотно-зависимыми характеристиками. Основой алгоритма является методика расчета собственных частот колебаний систем с вязкоупругими элементами, характеристики которых не зависят от частоты, предложенная в работе [35]. Она построена на использовании возможностей пакета ANSYS, позволяющих осуществлять запись глобальных ансамблированных матриц жесткости и масс, и метода Мюллера [33], который позволяет определять искомые комплексные собственные значения матриц, т.е. комплексные собственные частоты колебаний.

Для простоты и наглядности процесс построения алгоритма численного решения задачи о собственных колебаниях для конструкций, содержащих элементы из частотно-зависимых вязкоупругих материалов, будет рассмотрен на примере изотропного вязкоупругого материла, который характеризуется четырьмя компонентами комплексных динамических модулей: $G_{Re}$, $G_{Im}$, $B_{Re}$, $B_{Im}$.

Пусть для каждой составляющей комплексного динамического модуля известна ее графическая зависимость от частоты внешнего воздействия $\Omega$, аналогичная представленным на рис. 1, из которых видно, что каждому i-тому значению частоты внешнего воздействия ${\Omega }_i$ соответствует единственно возможный набор значений компонент комплексных динамических модулей $G_{Re}^{\left(i\right)}$, $G_{Im}^{\left(i\right)}$, $B_{Re}^{\left(i\right)}$,$B_{Im}^{\left(i\right)}$.

 

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

[Figure 1. Graphical relations between components of complex dynamic moduli and frequency of external excitation]

 

Зная значения компонент комплексных динамических модулей $G_{Re}^{\left(i\right)}$, $G_{Im}^{\left(i\right)}$, $B_{Re}^{\left(i\right)}$, $B_{Im}^{\left(i\right)}$, снятых с графика, можно определить значения комплексных собственных частот колебаний ${\left\{\omega \right\}}_i^{\top }={\left\{{\omega }_{Re}\right\}}_i^{\top }+i\{{{\omega }_{Im}\}}_i^{\top }$, соответствующих этим значениям материальных характеристик. Фигурные скобки в данном случае обозначают вектор-столбцы, содержащие значения комплексных собственных частот колебаний, соответствующих набору материальных констант $G_{Re}^{\left(i\right)}$, $G_{Im}^{\left(i\right)}$,$B_{Re}^{\left(i\right)}$, $B_{Im}^{\left(i\right)}$, полученных для i-того значения частоты внешнего воздействия ${\Omega }_i$. Таким образом, полученный вектор-столбец ${\left\{\omega \right\}}_i^{\top }={\left\{{\omega }_{Re}\right\}}_i^{\top }+i\{{{\omega }_{Im}\}}_i^{\top }$определяет спектр собственных частот колебаний, полученный для -того набора значений компонент комплексных динамических модулей.

Если вычислить ${\left\{\omega \right\}}_i^{\top }$для $\left(n\right)$ наборов материальных констант, то можно построить зависимость действительных ${\omega }_{Re}$ и мнимых ${\omega }_{Im}$ частей комплексных собственных частот колебаний от значений компонент комплексных динамических модулей. Обозначим эти зависимости как ${\omega }_{Re}\left(\Omega \right)$ и ${\omega }_{Im}\left(\Omega \right)$. При этом параметр $\Omega$, обозначающий частоту внешнего возбуждения, определяет соответствующий этой частоте набор материальных характеристик $G_{Re}$, $G_{Im}$, $B_{Re}$,$B_{Im}$. Таким образом, зависимости ${\omega }_{Re}\left(\Omega \right)$ и ${\omega }_{Im}\left(\Omega \right)$ есть зависимости ${\omega }_{Re}\left(G_{Re}\mathrm{,}{G}_{Im}\mathrm{,}{B}_{Re}\mathrm{,}{B}_{Im}\right)$ и ${\omega }_{Im}\left(G_{Re}\mathrm{,}{G}_{Im}\mathrm{,}{B}_{Re}\mathrm{,}{B}_{Im}\right)$.

В момент резонанса значение частоты внешнего возбуждения совпадает со значением собственной частоты колебаний объекта. При этом в случае, когда собственные частоты колебаний являются комплексными, резонанс наступает при совпадении со значением частоты внешнего возбуждения значения действительной части комплексной собственной частоты ${\omega }_{Re}$. Данное условие можно выразить следующим образом:

${\omega }_{Re}=\Omega . \left(17\right)$

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

 

Рис. 2. Зависимости действительных ω_Re (a) и мнимых частей ω_Im (b) комплексных собственных частот колебаний от значений компонент комплексных динамических модулей

[Figure 2. Relations between real parts ω_Re ( a) and imaginary parts ω_Im (b) of complex natural vibration frequency and values of component of complex dynamic moduli]

 

Рис. 3. Методика определения значений действительных ω_Re (a) и мнимых частей ω_Im (b) комплексных собственных частот колебаний

[Figure 3. A technique for determining values of real parts ω_Re (a) and imaginary parts ω_Im (b) of complex natural vibration frequencies]

 

После определения значений действительных частей комплексных собственных частот колебаний определяются соответствующие им значения мнимых частей. Для этого на график зависимости мнимых частей выбранных мод колебаний от частоты внешнего возбуждения (рис. 2, b) наносятся прямые, параллельные оси ординат, которые пересекают ось абсцисс в точках, соответствующих условию (17). Пересечение каждой из этих прямых с кривой мнимой части для соответствующей моды даст значение мнимой части комплексной собственной частоты колебаний (рис. 3, b).

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

4. Апробация и верификация алгоритма. Продемонстрируем предлагаемый алгоритм определения комплексных собственных частот колебаний вязкоупругой конструкции с частотно-зависимыми характеристиками материала на примере двухслойной консольно защемленной пластинки, представленной на рис. 4. Один слой пластинки с размерами $l_p=l_v=210$ мм, $b_p=b_v=26$ мм, $h_p=h_v=0.6$ мм выполнен из упругого материала (обозначен голубым цветом на рис. 4), а второй — из вязкоупругого (обозначен желтым цветом). Физико-механические характеристики упругого слоя следующие: модуль Юнга $E=2\cdot {10}^{11}$Па, коэффициент Пуассона $\nu =0.3$, удельная плотность $\rho =7800$ кг/м.

 

Рис. 4. Расчетная схема конструкции (онлайн в цвете)

[Figure 4 (color online). Computational scheme of the structure]

 

Пусть свойства вязкоупругого материала в диапазоне частот внешнего воздействия от 0 до 500 Гц изменяются так, как представлено на рис. 5. Данные зависимости качественно отражают один из возможных вариантов поведения компонент частотно-зависимых комплексных динамических модулей [21, 45] и выбраны исключительно для демонстрации работоспособности и эффективности предлагаемого алгоритма.

 

Рис. 5. Графические зависимости компонент комплексных динамических модулей от частоты внешнего воздействия: возрастающие (a) и убывающие (b)

[Figure 5. Graphical relations of components of complex dynamic moduli and frequency of external excitation: increasing (a) and decreasing (b)]

 

Диапазон изменения составляющих комплексных динамических модулей лежит в следующих пределах: $G_{Re}=1.5\cdot {10}^9÷1.26\cdot {10}^{10}$ Па, $G_{Im}=0.2G_{Re}$, $B_{Re}=7.45\cdot {10}^{10}÷6.26\cdot {10}^{11}$Па. Мнимая часть комплексного модуля объемного сжатия принята равной нулю $B_{Im}=0$. Удельная плотность материала вязкоупругого слоя равна ${\rho }_v=1200$кг/м.

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

Комплексные собственные частоты колебаний, полученные на основе решения задачи о собственных колебаниях, и резонансные частоты, определенные на амплитудно-частотных характеристиках (АЧХ) в результате решения задачи о вынужденных установившихся колебаниях в пакете прикладных программ ANSYS, сравниваются с полученными по предлагаемому алгоритму действительными составляющими собственных частот конструкции.

Рассмотрим вынужденные установившиеся колебания кусочно-однородных вязкоупругих тел, свойства которых описываются частотно-зависимыми комплексными динамическими модулями. При анализе вынужденных установившихся колебаний были рассмотрены оба варианта зависимости компонент комплексных динамических модулей от частоты (возрастающие с ростом частоты и убывающие), представленные на рис. 5. Возбуждение колебаний осуществлялось путем приложения к защемленному концу пластины вектора перемещений ${\overline{U}}_0=\left\{U_x\mathrm{,}{U}_y\mathrm{,}{U}_z\right\}=\{1,1,1\}$ мм.

На рис. 6 приведены первые 6 мод колебаний с указанием точки съема информации о компонентах вектора перемещений на данной моде (рис. 6, a); полученные амплитудно-частотные характеристики (АЧХ) модуля вектора перемещений $\left|U_{sum}\right|=\sqrt{U_x^2+U_y^2+U_z^2}$, отнесенных к модулю вектора возбуждающего усилия $\left|U_0\right|=\sqrt{U_x^2+U_y^2+U_z^2}$ для первых шести мод в рамках выбранного диапазона частот внешнего воздействия (рис. 6, b, c).

На рис. 6, b приведены результаты, полученные для компонент комплексных динамических модулей вязкоупругого материала, возрастающих с ростом частоты внешнего воздействия, а на рис. 6,  c — для убывающих. Значения компонент вектора перемещений снимались в точке  (рис. 6, a). Стоит отметить, что для гарантированной регистрации крутильной (№ 3, рис. 6, a) и планарной (№ 6, рис. 6, a) мод колебаний расположение точки A отличается от положения остальных точек съема АЧХ для изгибных мод колебаний. Это объясняется тем, что величина перемещений на изгибных модах колебаний существенно выше, чем на крутильных или планарных модах. В результате при фиксации изгибных и, например, крутильных мод в одной и той же точке при построении АЧХ резонанс крутильной моды может попадать на восходящую или нисходящую ветвь резонансной кривой изгибной моды и в силу большой разницы в величине перемещений, будет на графике не виден. Перенос же точки съема информации позволяет его обнаружить.

 

Рис. 6. Первые шесть мод колебаний пластинки (a) и амплитудно-частотные характеристики |U_sum|/|U₀|, полученные для возрастающих (b) или убывающих (c) с ростом частоты внешнего воздействия компонент комплексного динамического модуля

[Figure 6. The first six modes of vibrations of the plate ( a) and frequency response plots of |U_sum|/|U₀| obtained for increasing (b) and decreasing (c) components of complex dynamic moduli]

 

Стоит отметить, что в случае убывающих с ростом частоты значений комплексных модулей для рассматриваемой системы наблюдается смена порядка чередования 5 и 6 форм колебаний. Так, для убывающих модулей мода, соответствующая значению 5-й собственной частоты ${\omega }_5=423.29$ Гц, становится изгибной, а мода, соответствующая 6-й собственной частоте ${\omega }_6=450.84$ Гц, — планарной. Тогда как для возрастающих модулей мода, соответствующая значению 5-й собственной частоты ${\omega }_5=460.20$ Гц, является планарной, а мода, соответствующая 6-й собственной частоте ${\omega }_5=499.32$ Гц, — изгибной.

Пунктирной линией на рисунках 6, b и 6, c обозначены значения собственных частот колебаний, определенных по предложенному алгоритму.

В таблице приведено сравнение собственных частот колебаний конструкции, в составе которой есть слой, выполненный из материала с частотно-зависимыми возрастающими или убывающими составляющими комплексных динамических модулей. Для сравнения в столбцах 2 и 4 данной таблицы приведены значения частот колебаний, соответствующих резонансным пикам на АЧХ и полученных в результате решения задачи о вынужденных установившихся колебаниях и представленных на рис. 6.

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

 

Комплексные собственные частоты колебаний, определенные на основе предлагаемого алгоритма, и резонансные частоты, определенные на основе задачи о вынужденных колебаниях, решенной в ANSYS [Compex natural vibration frequencies detrmined on the basis of the proposed algorithm and resonant frequencies determined on the basis of forced steady-state vibration problem solved in ANSYS]

	Incremental moduli		Decreasing moduli
Modes of	ANSYS	Suggested	ANSYS	Suggested
vibrations		method		method
1	11.33	11.34 +  i 0.16	15.21	15.24 + i 0.72
2	77.11	77.09 + i 2.03	88.34	88.29 + i 3.63
3	203.66	206.26 + i 7.02	197.00	199.28 + i 5.97
4	239.95	239.99 + i 9.20	225.96	225.81 + i 7.19
5	460.26	460.20 + i 3.01	423.39	423.75 + i 11.14
6	499.32	499.36 + i 21.73	450.84	450.85 +  1.22

 

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

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

Сравнение значений действительных частей комплексных собственных частот, определенных из решения задачи о собственных колебаниях, демонстрирует хорошее согласование со значениями резонансных частот, определенных из решения задачи о вынужденных установившихся колебаниях в пакете прикладных программ ANSYS.

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

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

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

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

Конкурирующие интересы. Конфликты интересов отсутствуют. Номер лицензии ANSYS в «ИМСС УрО РАН» № 1064623.

Авторский вклад и ответственность. Все авторы принимали участие в разработке концепции статьи и в написании рукописи. Авторы несут полную ответственность за предоставление окончательной рукописи в печать. Окончательная версия рукописи была одобрена всеми авторами.

Финансирование. Работа выполнена в рамках государственного задания ПФИЦ УрО РАН по теме № АААА–А19–119012290100–8.

Competing interests. Authors have no competing interests. ANSYS license number no. 1064623, Institute of Continuous Media Mechanics UB RAS.

Authors’ contributions and responsibilities. Each author has participated in the article concept development and in the manuscript writing. The authors are absolutely responsible for submit the final manuscript in print. Each author has approved the final version of manuscript.

Funding. The work was supported with the State Assignment of Perm Federal Research Center of the Ural Branch of the Russian Academy of Sciences, no. AAAA–A19–119012290100–8.

About the authors

Dmitrii A. Oshmarin

Institute of Continuous Media Mechanics UB RAS

Email: oshmarin@icmm.ru
ORCID iD: 0000-0002-9898-4823
SPIN-code: 6084-5189
Scopus Author ID: 57041319000
ResearcherId: J-7906-2018
http://www.mathnet.ru/person122272

Juior Researcher; Dept. of Complex Problems of Deformable Solids Mechanics

1, Academician Korolev str., Perm, 614013

Natalya V. Sevodina

Institute of Continuous Media Mechanics UB RAS

Email: natsev@icmm.ru
ORCID iD: 0000-0001-9374-7135
SPIN-code: 1605-0002
Scopus Author ID: 15133373300
http://www.mathnet.ru/person73434

Cand. Techn. Sci.; Researcher; Dept. of Complex Problems of Deformable Solids Mechanics

1, Academician Korolev str., Perm, 614013

Nataliya A. Iurlova

Institute of Continuous Media Mechanics UB RAS

Author for correspondence.
Email: yurlova@icmm.ru
ORCID iD: 0000-0003-3497-0358
SPIN-code: 7391-6762
Scopus Author ID: 57191952953
ResearcherId: N-5129-2015
http://www.mathnet.ru/person122838

Cand Phys. & Math. Sci.; Associate Professor; Academic Secretary

1, Academician Korolev str., Perm, 614013

References

Supplementary files