Проблемы машиностроения и надежности машин, 2019, № 6, стр. 59-70

1. Решетчатые конструкции широко распространены в авиа- и судостроении, во многих других отраслях современного машиностроения, в строительстве и сельском хозяйстве. Важный подкласс таких систем – плоские 2D-решетки. Они определяют конструкции ряда горных, просеивающих и сортировочных машин, с их помощью проводится моделирование многочисленных 2D-объектов – мембран, пластин, панелей, строительных конструкций. Кроме того, при их посредстве могут анализироваться волны в метаматериалах, композитных материалах с периодическим армированием, кристаллах, а также наноструктурированных поверхностных и приповерхностных слоях конструкционных материалов [1–7]. Этот перечень можно существенно продолжить. Вместе с тем, динамический анализ решетчатых конструкций по ряду причин не получил необходимого распространения, а работы, касающиеся динамики решетчатых конструкций с учетом соударений их узлов достаточно редки [8–10].

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

2. Рассмотрим струнную решетку, составленную из трех семейств упругих линейных струн, образующих правильные треугольные ячейки (рис. 1a). Выберем на решетке декартову остроугольную систему координат, образованную осями x и y, расположенными друг к другу под углом π/3 (рис. 1б). Для идентификации узлов решетки, очевидно, достаточно двух текущих индексов. Обозначим их k = 0, 1, 2, …, N₁ и q = 0, 1, 2, …, N₂, предположив что индекс k нумерует вдоль оси x, а индекс q – оси y.

Пусть узлы – суть абсолютно твердые тела с массами ${{m}_{{kq}}} \equiv m$. Будем искать их перемещения (прогибы) вдоль осей, расположенных перпендикулярно плоскости решетки.

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

Обозначим натяжения каждого пролета струны τ, а длины сторон треугольников – a. Тогда коэффициент упругости с = τa^–1.

Каждый узел, находящийся внутри решетки имеет шесть равнозначных соседних узлов, поэтому систему можно трактовать как гексагональную центрированную решетку [11, 12]. Функция Лагранжа решетки имеет вид

где ${{T}_{{\text{*}}}}{\text{(}}{{\dot {u}}_{{kq}}}{\text{)}}$ и ${{U}_{{\text{*}}}}{\text{(}}{{u}_{{kq}}}{\text{)}}$ – кинетическая и потенциальная энергии; $D{}_{n}^{{}}{\text{(}}u_{{kq}}^{{}}{\text{)}}$ – изменения длин струнных пролетов, между узлом (k, q) и его шестью соседями (n = 1, …, 6).

Пусть в узлах приложены внешние силы, представимые через операторные функции g_kq($d{\text{/}}dt$; t, $u_{{kq}}^{{}}$). Тогда для изгибных колебаний решетки получаем уравнения движения

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

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

Рассматриваемая система аналогична решетчатой конструкции с квадратными (прямоугольными) ячейками (рис. 2а) [8–10]. Уравнения движения для прогибов решетки с квадратными ячейками имеют вид

где a – длина стороны квадратной ячейки. Вид граничных условий (3) – сохраняется.

3. Перейдем к континуальному аналогу задачи (2), (3). Предположим, что число струн, образующих решетку велико, т.е. ее локальная структурная характеристика (длина стороны треугольной ячейки a) малая величина. Одновременно усредненные характеристики решетки (поверхностная плотность (масса единицы площади) $\rho = m{\text{/}}{{a}^{2}}$ и упругость пролетов с) не малы. Постулируем предположение, что существует требуемое число раз непрерывно-дифференцируемая функция $u(x,y,t)$, совпадающая в узлах решетки с дискретной физической переменной ${{u}_{{kq}}} = u{\text{(}}x,y,t{\text{)}}$. При этом, с точностью до бесконечно малых ($O{\text{(}}{{a}^{2}}{\text{)}}$)

В то же время, можно переобозначить ${{g}_{{kq}}}{\text{(}}p;t,{{u}_{{kq}}}{\text{)}}$ = $g{\text{[}}x,y;p;t,u{\text{(}}x,y,t{\text{)]}}$. Подставив (5) в (2), проведем очевидные разложения в степенные ряды и получим: где ${{g}_{{\text{0}}}} = g{{a}^{{ - 2}}}$.

Граничные условия (3) переходят в условия

где P контур, на котором находятся, защемленные граничные узлы решетки – параллелограмм или какой-то иной контур защемления.

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

Предположим, что решетка не ограничена и ${{g}_{0}} = 0$. Представим решения уравнений (6) или (8) в виде синусоидальных волновых пакетов где {K₁, K₂} – компоненты волнового вектора; ω – частота. Отсюда можно получить дисперсионные соотношения: для системы (8)

Для системы (6)

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

Отметим, что длинноволновые приближения для уравнений движения решетки с квадратными ячейками (4) дадут вместо (8) и (6) соответственно уравнения мембраны

При этом дисперсионные соотношения принимают вид

Для основной дискретной модели треугольной решетки, полагая

после ряда преобразований, следующих из внесения (15) в уравнение (2), найдем дисперсионное соотношение вида Для длинных волн ${{K}_{1}}$, ${{K}_{2}} \to 0$ и из (16) находим что с точностью до множителя совпадает с главной частью дисперсионного соотношения (11).

Рассмотрим вопрос о том, как соотносятся между собой модели решеток с треугольными и квадратными ячейками. Возьмем модели, отвечающие уравнениям (8) и (12). Их можно объединить посредством введения некоторого безразмерного множителя перед смешанной производной (обозначим его $\sigma $). В квадратной скобке, входящей в левую часть уравнения (8) появляется член $\sigma u_{{xy}}^{{}}(x,y,t)$. Пусть $\sigma \in {\text{[0,1]}}$. Если $\sigma $ = 1, то имеем модель (8), при $\sigma \to 0$ приходим к модели (12).

Если ввести тот же параметр $\sigma $ перед разностями ${\text{(}}{{u}_{{kq}}} - {{u}_{{(k + 1,q + 1)}}}{\text{)}}$ и (${{u}_{{kq}}} - {{u}_{{(k - 1,q - 1)}}})$, то модель решетки с треугольными ячейками (уравнение (2)) перейдет в модель с квадратными ячейками (уравнение (4)).

Корректное описание преобразования решетки при изменении параметра $\sigma $ – нетривиально, так как несмотря на то, что треугольная и прямоугольная ячейки – топологически эквивалентны [13] (рис. 3) и могут быть непрерывно деформированы одна в другую, придать однозначный физический смысл промежуточным значениям параметра $\sigma \in {\text{(0,1)}}$ без постулирования дополнительных гипотез проблематично. Относящиеся сюда вопросы – не рассматриваются.

4. Дальнейшие рассмотрения касаются только модели (2) и случая вынужденных случайных колебаний, происходящих в вязкой среде. Будем предполагать, что вблизи каждого из узлов решетки с зазором $\Delta \geqslant 0$ установлен, вообще говоря, податливый ограничитель хода (рис. 4; показана одна гексагональная ячейка решетки). Пусть возбуждение колебаний осуществляется за счет действия на все узлы случайных сил, описываемых стационарными стандартными белыми шумами с нулевым средним. Взаимодействие узлов решетки с ограничителями предполагается упругим. Виброударный вариант модели (4) (рис. 2б) рассматривался ранее в работах [8–10].

Будем использовать канонические переменные ${{u}_{{kq}}} \equiv {{q}_{{kq}}}$; ${{p}_{{kq}}} = m{{\dot {u}}_{{kq}}}$. Пусть $\lambda {\text{П(}}q_{{kq}}^{{}}{\text{)}}$ потенциальная энергия силы упругого ударного взаимодействия, определяемого гипотезой удара, $\lambda \gg 1$ – большой параметр. Считая, что принятые гипотезы взаимодействия (удара) во всех ударных парах одинаковы, сможем записать

где $\lambda {\text{Ф(}}q_{{kq}}^{{}}{\text{)}}$ – сила взаимодействия.

Установление гипотезы удара – это мотивированное задание представления силы взаимодействия. Если, например, ударная пара – односторонняя и несимметричная, то несимметричной будет и, отвечающая ей, пороговая функция $\lambda {{{\text{Ф}}}_{{kq}}}({{q}_{k}})$ (рис. 4), при посредстве которой, собственно, и устанавливается гипотеза удара в данной ударной паре. Класс пороговых функций {Ф}_Δ для данной задачи определяется так

Здесь предполагаем, что функция ${\alpha (}q{}_{{kq}}{\text{)}}$ – непрерывно дифференцируема на всей числовой оси, монотонно возрастает и является выпуклой при ${{q}_{{kq}}}$ ≥ 0 (рис. 4); η(х) – единичная функция. Если взаимодействие – двустороннее, то функция взаимодействия будет двухпороговой. После первоначальной работы [14] подобный подход получил развитие в книге [15], а для виброударных систем со многими ударными парами, в частности, для 2D-систем – в работах [9, 15 и др.].

При $\lambda \to \infty $ гипотеза взаимодействия переходит в гипотезу Ньютона для абсолютно упругого удара. Т.е. конфигурационное пространство системы $\Re $ задается при помощи неравенств

При соударениях в какие-то моменты времени ${{t}_{\alpha }}:{{p}_{{kq}}}({{t}_{\alpha }} - 0)$ = $ - {{p}_{{kq}}}({{t}_{\alpha }} + 0)$.

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

Здесь по условию: ${{\xi }_{{kq}}}{\text{(}}t{\text{)}}$ – случайные силы, описываемые стандартными белыми шумами с одинаковыми интенсивностями 2S: $\left\langle {{{\xi }_{{kq}}}{{\xi }_{{rs}}}} \right\rangle $ = $2S{{\delta }_{{kq}}}{{\delta }_{{rs}}}\delta (t - t{\text{')}}$; ${{\delta }_{{kq}}}$, ${{\delta }_{{rs}}}$ – символы Кронекера (k, r =1, 2, … N₁; q, s = 1, 2, …, N₂) угловые скобки обозначают операцию статистического усреднения; δ(t) – δ – функция.

Для консервативной части системы функцию Гамильтона $H({{q}_{{kq}}},{{p}_{{kq}}})$ легко определить, исходя из вида функции Лагранжа (1) и формулы (17), а также предположения, что потери энергии при соударениях в каждой ударной паре не учитываются

Рассматриваемая решетка при любых значениях чисел ${{N}_{1}}$, ${{N}_{2}}$, определяющих ее размеры будет отвечать уравнениям Гамильтона вида

Воспользуемся методами диффузионных марковских процессов и запишем уравнение ФПК для w(t, q_kq, p_kq) – совместной плотности вероятностей многомерного случайного процесса {q_kq(t), p_kq(t)} [16–19], которая и определяет состояние решетки

– скобка Пуассона для узла (k, q); A и B – произвольные гладкие функции. Каждая компонента матриц q_kq(t) и p_kq(t) дает соответствующее значение канонических переменных q_kq и p_kq.

Как известно [16], для стационарного процесса, когда $\frac{{\partial P}}{{\partial t}} = 0$, уравнение (6) имеет точное решение, которое дается каноническим распределением Гиббса

в чем можно убедиться, непосредственно, подставив формулу (22) в (21) и воспользоваться условием стационарности; константа C определяется после нормировки.

Хорошо известно (и это видно из формул (19) и (22)), что импульсы и координаты – статистически независимы.

Распределение Гиббса после интегрирования по всем координатам дает нормальное распределение по импульсным переменным (распределение Максвелла): кинетическая энергия, стоящая в показателе экспоненты, зависит только от квадратов импульсов

Для того, чтобы установить его нормальность достаточно перенумеровать импульсы и перейти от двойной к последовательной индексации. С помощью формулы (23) можно получить всю необходимую информацию об импульсах и силовых воздействиях в системе [15, 19].

5. Обратимся, после интегрирования распределения (22) по всем обобщенным импульсам, к распределению Больцмана по обобщенным координатам

где учтен вид функции Гамильтона (19).

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

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

Функция распределения координат

где подразумевается, что проводится N₁N₂-кратное интегрирование по всем координатам и символ $d{{{\mathbf{z}}}_{{kq}}}$ – суть произведение N₁N₂ – дифференциалов; ${\mathbf{P}}{\text{\{ }}{{{\mathbf{q}}}_{{\text{*}}}}_{{kq}} < {{{\mathbf{q}}}_{{kq}}}{\text{\} }}$ – вероятность события, при котором покомпонентно имеет место указанное отношение порядка (<).

Отсюда видно, что вероятность существования конфигурации с превышением порога $\Delta $ для всех узлов решетки есть

где матрица Δ отвечает достижению всеми узлами решетки пороговой конфигурации.

Внося (25) в (26), и далее, внося результат в (29) находим

При больших λ оценка вероятности (28) есть [21]

В то же время, можно показать, что для любых фиксированных значений индексов a и b вероятность пребывания за порогом единственного узла решетки P{u_ab > $\Delta $} = = O(${{{\lambda }}^{{ - {\text{ }}\frac{1}{2}}}}$), для пар, троек, “энок” при любых допустимых значениях индексов при $\lambda \to \infty $ получается: P{u_ab > $\Delta $; u_cd > $\Delta $} = O(${{\lambda }^{{ - 1}}}$); P{u_ab > $\Delta $; u_cd > $\Delta $; u_cd >} = O(${{\lambda }^{{ - \frac{3}{2}}}}$); P{u₁ > > $\Delta $; …; u_n > $\Delta $} = O(${{\lambda }^{{ - \frac{n}{2}}}}$). Таким образом реализация конфигураций, сопровождающихся пребыванием точечных тел за пороговыми значениями маловероятны. Эти асимптотики могут быть уточнены и переписаны в размерных формах при помощи общих приемов методов асимптотических оценок [21].

Если ${\lambda } \to \infty $, то высокие при конечных $\lambda $ потенциальные барьеры, $\lambda {\text{П(}}q_{{kq}}^{{}}{\text{)}}$ приобретают высоты равные кинетическим энергиям каждого из узлов и становятся непроницаемыми [9, 14]. Времена взаимодействий обращаются в пределе в ноль, а гипотеза удара переходит в гипотезу Ньютона для упругого взаимодействия.

При этом вид распределения Гиббса изменится. В представлении (23) исчезает член с большим параметром в показателе экспоненты и добавляются соотношения на границе и условия удара. В данном случае

причем два последних соотношениях имеют место для любого узла (k, q); t_0kq – произвольный момент его соударения.

Распределение Максвелла (23) для случайных обобщенных импульсов сохраняет свой вид и удовлетворяет приведенным условиям упругого удара. Для анализа важны различного рода характеристики случайных импульсов ударов J_kq = 2 | p_kq |. После вычислений, для многомерного случайного процесса пользуясь работами [8, 22] приходим к усеченному нормальному распределению для случайного процесса J_k

Распределение (31) определяет силовые эффекты, связанные с влиянием ударов на динамику многомерной конструкции.

Для нечетных и четных j-х моментов одномерных величин J_kq можно, получить формулы [15, 19]

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

Рассматривая случайные конфигурации, возникающие при реализации виброударных режимов в решетке, примем во внимание малую вероятность “запороговых” конфигураций (см. (29)). При $\lambda \to \infty $ из (24) и (30) найдем, после интегрирования по импульсам, распределение Больцмана

где смысл обозначений станет очевиден из соотношений (1) и (19) и предполагается, что индексы k и q принимают все допустимые значения.

Получив представления (30) и (32), можно полностью описать “конфигурационные” характеристики системы: типы случайных профилей, их статистические характеристики и др.

Предположим, что есть возможность каким-то образом регистрировать любой случайно реализуемый профиль ${\mathbf{q}} = {\text{\{ }}q_{{kq}}^{{}}{\text{\} }}$. Зафиксируем его ${\mathbf{q}} = {{{\mathbf{q}}}^{{\text{0}}}}$ ≡ $\{ q_{{kq}}^{0}\} $ и определим среднюю частоту $({{\Omega }_{0}}\{ q_{{kq}}^{0}\} )$ его появления как математическое ожидание пересечения узлами решетки соответствующей поверхности в конфигурационном пространстве с положительной скоростью. Так как конфигурация $\{ q_{{kq}}^{0}\} $ – фиксирована, то фиксировано и значение потенциальной энергии ${{U}_{{\text{*}}}}(q_{{kq}}^{{\text{0}}})$ и, следовательно, можно найти

Легко увидеть, что при функционировании рассматриваемой системы наиболее часто будет регистрироваться равновесная конфигурация ${{{\mathbf{q}}}^{{\text{0}}}} \equiv {\mathbf{0}}$.

Так как ${{U}_{{\text{*}}}}{\text{(}}{\mathbf{0}}{\text{)}} = {\text{0}}$, то после интегрирования из формулы (33) найдем

Рассмотрим случайные профили решетки, при появлении которых в некоторые моменты времени, узлы либо находятся в равновесном состоянии, либо выброшены на ограничитель. Обозначим произвольную конфигурацию ${\mathbf{q}}{\text{(}}\Delta {\text{;}}\beta {\text{)}}$ = ${\text{\{ }}q_{{kq}}^{{}}{\text{(}}\Delta {\text{;}}\beta {\text{)\} }}$, где индекс $\beta \in {\mathbf{B}}$ нумерует описанные профили; число элементов множества B, очевидно, равно ${{2}^{{{{N}_{1}}{{N}_{2}}}}}$.

Из формул (33) и (34) найдем

Обозначим среднюю относительную частоту появления рассматриваемой конфигурации ${{\omega }_{\beta }}$ = ${{\omega }_{\beta }}{\text{\{ }}q_{{kq}}^{{}}{\text{(}}\Delta {\text{;}}\beta {\text{)\} }}$ = ${{\Omega }_{0}}{\text{\{ }}q_{{kq}}^{{}}{\text{(}}\Delta {\text{;}}\beta {\text{)\} /}}{{\Omega }_{0}}{\text{\{ }}{\mathbf{0}}{\text{\} }}$. Из формулы (35) следует

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

Рассматриваемые профили отвечают двум характерным значениям обобщенных координат: нахождение в положении равновесия и достижение ограничителей: ${\mathbf{q}}{\text{(}}\Delta {\text{;}}\beta {\text{)}}$ = = . Они, очевидно, достаточно характерны для данных типов решеток и их удобно регистрировать при помощи диаграмм показанных на рис. 5, где выбрана решетка в форме параллелограмма, имеющего 28 узлов, включая лежащие на границе. Будем выделять узлы, находящиеся у ограничителя как показано на рис. 5а. В случае пребывания в состоянии равновесия выделять узел не будем. Таким образом можно получить диаграммы, отвечающие рассматриваемым конфигурациям при помощи которых можно оценить величины соответствующих значений потенциальной энергии, фигурирующей в формуле (36).

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

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

Видно, что чем выше “кучность” соударяющихся узлов, тем меньше будет соответствующий уровень потенциальной энергии. Из диаграмм видно, что упругие элементы соседних взаимодействующих узлов не деформируются. Таким образом соседство двух выделенных или, наоборот, не выделенных узлов не увеличивает значение потенциальной энергии конфигурации. В то же время соседство выделенного и не выделенного узла увеличивает потенциальную энергию на величину ~$с{{\Delta }^{2}}$. Так что для конфигурации диаграммы b имеем 16 деформаций и ${{U}_{{\text{*}}}}$ ~ 16$с{{\Delta }^{2}}$, для диаграммы c – ${{U}_{{\text{*}}}}$ ~ ~ 30$с{{\Delta }^{2}}$, а для некоторой промежуточной конфигурации d – ${{U}_{{\text{*}}}}$ ~ 20$с{{\Delta }^{2}}$. Получается, что относительная частота появления конфигурации b – наиболее высока для комбинации с шестью соударяющимися узлами.

7. Приведем заключительные замечания.

7.1. Отмечено, что наименьшие значения потенциальной энергии имеют конфигурации, состоящие из единственного соударяющегося узла. Простые подсчеты показывают, что при рассмотрении профилей, в которых соударяются $Q \geqslant 2$ узлов, наименьшей потенциальной энергией обладают те, которые составляют выпуклый многоугольник. При $Q = 2$ многоугольник вырождается в отрезок.

7.2. Вместо неравенств (18) в рамках данной модели допустимо рассмотрение более сложной структуры конфигурационного пространства

Это означает, что каждый узел может быть снабжен собственным ограничителем. При этом бесконечные значения отвечают случаю его отсутствия для данного узла (k, q). Распределения (30) и (31) – сохраняются, а для распределения (32) фазовые ограничения видоизменятся в соответствии с неравенствами (37). Структура формулы (36) сохранится.

7.3. При рассмотрении решеток с квадратными (прямоугольными) ячейками [8] все качественные выводы сохраняются. При рассмотрении одномерных цепочек [22] относительно часто регистрируемыми оказываются профили с одиночными, двойными, тройными и т.д. выбросами на ограничитель. При этом либо узел единственный, либо все узлы – лежат на одном отрезке. В этом – принципиальное отличие одномерной модели.

Рассмотренные модели могут быть использованы при рассмотрении двумерных дискретно-континуальных систем [22], а также при изучении самоподобных решетчатых структур [23].

7.4. Инженерное значение рассмотренных задач отмечено в п. 1. Они имеют не только модельную ценность, но и представляют интерес для инженеров-расчетчиков, так как открывают возможность получить определенные предварительные представления при выполнении динамических расчетов ряда конкретных двумерных конструкций, функционирующих в условиях соударения их элементов. В данном случае задачи усложняются и необходимо привлечение численных методов [24].

Информация о финансовой поддержке. Работа выполнена за счет гранта РНФ (проект № 19-19-00065).