Spatial Dispersion of Acoustic Waves in Functionally Graded Rods

1. Введение. 1.1. Акустические волны в ФГ-средах. Последние исследования поверхностных акустических волн (ПАВ) в функционально-градиентных (ФГ) и слоистых пластинах с поперечной неоднородностью выявили ряд особенностей дисперсионных свойств, так как точки с нулевой групповой скорости (НГВ), принадлежащие фундаментальным ветвям, точки пересечения и расхождения между фундаментальной ветвью и нижними ветвями высших мод и др. [1-10]. Большое количество работ посвящено развитию численных методов решения дисперсионных уравнений для направленных волн в ФГ-пластинах; так, в [11] используются ряды Пеано; проекционный метод Галеркина предложен в [12]; несколько подходов связаны с разложением упругих потенциалов, полей напряжений и перемещений и варьированием свойств материала по полиномам Лежандра [13-15] или применением метода спектральных элементов (СЭМ); см. [16, 17].

Дисперсионные волны Похгаммера$–$Кри в слоистых и ФГ-стержнях с поперечной неоднородностью изучаются либо с помощью аналитических подходов [18$–$22], подходящих для некоторых конкретных видов неоднородностей, либо численно, чаще всего с помощью СЭМ в пользу МКЭ; см. [23] или полиномы Лежандра [24]. Хотя волны в одномерных ФГ-стержнях могут не иметь таких обильных дисперсионных явлений, которые наблюдаются в ФГ-пластинах [25] и трехмерных или аксиально-симметричных ФГ-стержнях, тем не менее, как показано ниже, волны в одномерных стержнях обладают несколькими необычными свойствами в отношении дисперсии и изменения энергии. Достаточно большое количество работ посвящено волновой динамике ФГ-сред с пьезоэлектрическими [26, 27], пористыми пьезоэлектрическими [28, 29] и пьезомагнитными свойствами [30].

1.2. Предмет исследования. Настоящее исследование направлено на: (i) построение определяющего уравнения и его решение для акустических волн в полубесконечном ФГ-стержне 1D с произвольной продольной неоднородностью; (ii) определение пространственного изменения энергии и величины осциллирующих волн; и (iii) получение дисперсионных свойств осциллирующих волн. Предполагается, что пространственная неоднородность относится к классу С¹(R₊), т. е. непрерывно дифференцируема в интервале x ∈ (0;+ ∞); см. рис. 1. Предполагается также, что соответствующие физические свойства удовлетворяют условию сильной эллиптичности по $R_{+}$ и имеют конечный и ненулевой пределы при $x\to +0$ и x → ∞; более подробно эти условия обсуждаются в следующем разделе. Решение в замкнутой форме для дисперсионного уравнения строится комбинацией шестимерного формализма Коши и экспоненциального метода фундаментальных матриц, оба метода ранее были разработаны для дисперсионного анализа направленных волн в слоистых и ФГ-пластинах [25].

Рис. 1. 1D ФГ-стержень; n показывает направление нормали волны и x показывает направление координатной оси вдоль направления распространения волны.

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

2. Основные уравнения. 2.1. Уравнение движения. Уравнение движения одномерного упругого стержня с продольной неоднородностью класса $C^1(R_{+})$ имеет вид [31]

         $\left(\rho (x){\partial }_{tt}^2-{\partial }_{x}E(x){\partial }_x\right)u(x,t)=0$,                   (2.1)

где $\rho (x)$ $–$ плотность материала; $E(x)$ $–$ модуль упругости; $u(x,t)$ $–$ продольное перемещение; $x$ $–$ пространственная переменная; $t$ $–$ время. Для обеспечения гиперболичности уравнения (2.1) накладываются следующие естественные ограничения

         $\rho (x)>\mathrm{0,}E(x)>0$,                      (2.2)

Чтобы представить уравнение в безразмерной форме, введем безразмерные постоянные ${\rho }_d,E_d,x_d,t_d$, имеющие размерность плотности, напряжений, длины и времени соответственно, дающие соответствующие безразмерные функции

${\rho }^{\prime }\left(x^{\prime }\right)=\frac{\rho \left(x\right)}{{\rho }_d}, E^{\prime }\left(x^{\prime }\right)=\frac{E\left(x\right)}{E_d}, u^{\prime }\left(x^{\prime }, t\text{'}\right)=\frac{u\left(x,t\right)}{x_d}$             (2.3)

где

         $x^{\prime }=\frac{x}{x_d}, t^{\prime }=\frac{t}{t_d}$.                              (2.4)

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

Введение новой переменной (напряжение)

         $\sigma (x,t)=E(x){\partial }_{x}u(x,t)$.                       (2.5)

Уравнение (2.1) превращается в систему дифференциальных уравнений первого порядка по переменной $x$

$\begin{array}{l}{\partial }_{x}u(x,t)=E{\left(x\right)}^{-1}\sigma (x,t)\\ {\partial }_x\sigma (x,t)=\rho \left(x\right){\partial }_{tt}^{2}u(x,t)\end{array}$.                   (2.6)

Следуя [32], рассмотрим гармоническую во времени волну с пространственно изменяющимися амплитудами

         $u\left(x,t\right)=m\left(x\right)e^{i\omega t},\sigma \left(x,t\right)=\tau \left(x\right)e^{i\omega t}$,                   (2.7)

где $m(x)$ $–$ искомая амплитуда, а $\omega$ $–$ безразмерная круговая частота.

Вводя новую векторную функцию

         $a(x)=\left(\begin{array}{c}m(x)\\ \tau (x)\end{array}\right)$                         (2.8)

и подстановка выражения (2.7) в уравнение (2.6) дает следующее матричное дифференциальное уравнение первого порядка

         ${\partial }_{x}a(x)=G(x)\cdot a(x)$,                         (2.9)

где

$G(x)=\left(\begin{array}{cc}0& E{(x)}^{-1}\\ -{\omega }^2\rho (x)& 0\end{array}\right)$.                   (2.10)

Уравнение (2.9) будем использовать для построения общих решений.

Рассматриваемые акустические волны, распространяющиеся в полубесконечном стержне, удовлетворяют следующим начальным и граничным условиям

${\left.u\left(x,t\right)\right|}_{t=+0}=\mathrm{0,} {\left.{\partial }_{t}u\left(x,t\right)\right|}_{t=+0}=\mathrm{0,} {\left.u\left(x,t\right)\right|}_{t=+0}=u_0{e}^{i\omega t}$.             (2.11)

Отметим, что с учетом гармонического во времени колебания (2.7) условие затухания Зоммерфельда становится менее ограничительным, чем для негармонических колебаний [33]:

         ${\left|u(x,t)\right|}_{x\to +\infty }<\infty$.                         (2.12)

Ввиду уравнения (2.7) обеспечение периодичности поля перемещений по переменной во времени, не существует предельных значений для функции $u(x,t)$ при $t\to \infty$ и какого-либо конечного x.

2.2. Общее решение. Общее решение уравнения (2.9) допускает следующий вид [34, 35]

         $a(x)=\exp \left(F(x)\right)\cdot a_0$,                         (2.13)

где $F(x)$ является первообразной $G(x)$ и $a_0$ $–$ вектор, определяемый граничными условиями.

         $F(x)={\displaystyle \int G(x^{\prime })d{x}^{\prime }}$.                           (2.14)

Предполагая матрицу $F$ полупростой и учитывая ее жорданову нормальную форму, получаем

         $F(x)=W(x)\cdot D(x)\cdot W^{-1}(x)$,                       (2.15)

где $W$ $–$ матрица, составленная из собственных векторов матрицы $F$, составленная по столбцам, и $D$ $–$ диагональная матрица, составленная из собственных векторов матрицы $F$. При разложении (2.15) матрица показателей в (2.13) имеет вид

         $\exp \left(F(x)\right)=W(x)\cdot \exp \left(D(x)\right)\cdot W^{-1}(x)$.                 (2.16)

С учетом уравнений (2.10), (2.14), матрицы $W(x)$ и $D(x)$ в представляются в следующем виде

         $W\left(x\right)=\left(\begin{array}{cc}\frac{\sqrt{a\left(x\right)}}{\sqrt{b\left(x\right)}}& -\frac{\sqrt{a\left(x\right)}}{\sqrt{b\left(x\right)}}\\ 1& 1\end{array}\right), D\left(x\right)=\left(\begin{array}{cc}\sqrt{a\left(x\right)b\left(x\right)}& 0\\ 1& -\sqrt{a\left(x\right)b\left(x\right)}\end{array}\right)$,       (2.17)

где

           $a\left(x\right)=\int \frac{dx\text{'}}{E\left(x^{\prime }\right)}, b\left(x\right)=-{\omega }^2\int \rho \left(x^{\prime }\right)dx\text{'}$.                  (2.18)

Отметим, что собственные векторы, расположенные в столбцах матрицы W, не нормализованы. Случай неполупростого $F$ приводит к появлению жордановой клетки в разложении и одного собственного и одного обобщенного собственных векторов в матрице W. Хотя такой случай физически возможен, он очень редок и в данной работе не рассматривается.

Выполняя операцию умножения в разложении (2.16) и учитывая (2.17)

         $\exp \left(D(x)\right)=\left(\begin{array}{cc}\exp \left(\sqrt{a(x)b(x)}\right)& 0\\ 0& \exp \left(-\sqrt{a(x)b(x)}\right)\end{array}\right)$           (2.19)

получаем

         $\exp \left(F(x)\right)=\left(\begin{array}{cc}\cosh \left(\sqrt{a(x)b(x)}\right)& \frac{\sqrt{a(x)}}{\sqrt{b(x)}}\sinh \left(\sqrt{a(x)b(x)}\right)\\ \frac{\sqrt{b(x)}}{\sqrt{a(x)}}\sinh \left(\sqrt{a(x)b(x)}\right)& \cosh \left(\sqrt{a(x)b(x)}\right)\end{array}\right)$.       (2.20)

Матрица показателей в (2.20) обладает интересным свойством.

Предложение 2.1. При любых E(x), ρ(x), удовлетворяющих условиям (2.2), матрица показателей $\exp \left(F(x)\right)$ не является вырожденной, более того

         $\det \left(\exp \left(F(x)\right)\right)=1$.                        (2.21)

Доказательство вытекает из рассмотрения правой части уравнения (2.20).

Следствие [36]. Ввиду уравнений (2.10), (2.14) и условия (2.21) экспоненциальное отображение $\exp \left(F(x)\right)$ сохраняет площадь.

Отметим, что ввиду отрицательного значения $b(x)$ в уравнении компоненты матрицы (2.17) являются мнимыми. Таким образом, матрица показателей в может быть записана через тригонометрические функции

         $\exp \left(F(x)\right)=\left(\begin{array}{cc}\cos \left(\sqrt{a(x)\left|b(x)\right|}\right)& \frac{\sqrt{a(x)}}{\sqrt{\left|b(x)\right|}}\sin \left(\sqrt{a(x)\left|b(x)\right|}\right)\\ -\frac{\sqrt{\left|b(x)\right|}}{\sqrt{a(x)}}\sin \left(\sqrt{a(x)\left|b(x)\right|}\right)& \cos \left(\sqrt{a(x)\left|b(x)\right|}\right)\end{array}\right)$.      (2.22)

2.3. Поля напряжений и перемещений. Граничные условия позволяют записать вектор $a_0$ из (2.13) следующим образом

         $a_0=\left(\begin{array}{c}{u}_0\\ 0\end{array}\right)$.                        (2.23)

Соединяя (2.13), (2.22) и (2.23), запишем выражения для полей напряжений и перемещений

$u\left(x,t\right)=cos\left(\left(\sqrt{a\left(x\right)\left|b\left(x\right)\right|}\right)\right)u_0{e}^{i\omega t}, \sigma \left(x,t\right)=-\sqrt{\frac{\left|b\left(x\right)\right|}{a\left(x\right)}}sin\left(\left(\sqrt{a\left(x\right)\left|b\left(x\right)\right|}\right)\right)u_0{e}^{i\omega t}$.      (2.24)

Откуда получим величины напряжений и перемещений

$\left|m\left(x\right)\right|=\left|cos\left(\sqrt{a\left(x\right)\left|b\left(x\right)\right|}\right)u_0\right|$, $\left|\tau \left(x\right)\right|=\left|\sqrt{\frac{\left|b\left(x\right)\right|}{a\left(x\right)}}sin\left(\left(\sqrt{a\left(x\right)\left|b\left(x\right)\right|}\right)\right)u_0\right|$.     (2.25)

2.4. Кинетическая и потенциальная энергия. Введем безразмерный период

         $T=\frac{2\pi }{\omega }$.                              (2.26)

Выражения (2.24) дают удельную кинетическую энергию и энергию деформации за период T:

         $\begin{array}{l}{E}_k\equiv \frac{1}{2}\rho (x){\displaystyle \underset{0}{\overset{T}{\int }}{\left|{\partial }_{t}u(x,t)\right|}^{2}dt}=\frac{\pi }{\omega }\rho (x){\cos }^2\left(\sqrt{a(x)\left|b(x)\right|}\right)u_0^2\\ E_s\equiv \frac{1}{2}{\left(E(x)\right)}^{-1}{\displaystyle \underset{0}{\overset{T}{\int }}{\left|\sigma (x,t)\right|}^{2}dt}=\frac{\pi }{\omega E(x)}\frac{\left|b(x)\right|}{a(x)}{\sin }^2\left(\sqrt{a(x)\left|b(x)\right|}\right)u_0^2\end{array}$.         (2.27)

3. Полиномиальные неоднородности. 3.1. Общие биномы. Класс положительных полиномов R может быть представлен в виде суммы двух квадратов [37].

$\rho \left(x\right)={\left({\displaystyle \sum }_{k=1}^n{\rho }_k{x}^k\right)}^2+{\rho }_0, E\left(x\right)={\left({\displaystyle \sum }_{k=1}^n{E}_k{x}^k\right)}^2+E_0$               (3.1)

с реальными коэффициентами ${\rho }_k,k=\mathrm{1,...,}n;E_k,k=\mathrm{1,...,}m$ и положительными ${\rho }_0$ и $E_0$.

Здесь рассматривается менее ограничительный класс положительных многочленов от R₊, это генерируется следующими биномами [38]

         $\rho \left(x\right)={\rho }_n{x}^n+{\rho }_0, E\left(x\right)=E_m{x}^m+E_0$                      (3.2)

с неотрицательными коэффициентами ${\rho }_n,E_m\ge 0$ и положительными ${\rho }_0,E_0>0$.

Подставляя биномы (3.2) в уравнение (2.18), получаем

$a\left(x\right)=\frac{1}{m{E}_0}Ф\left(-\frac{E_m}{E_0}{x}^m;1;\frac{1}{m}\right), b\left(x\right)=-{\omega }^2\left(\frac{{\rho }_n}{n+1}{x}^{n+1}+{\rho }_{0}x\right)$.        (3.3)

где $\Phi$ - $\Phi$ -функция Лерча, см. [39, 40]. Учитывая эти выражения, амплитуды в (2.25) будут выглядеть так:

$\begin{array}{l}\left|m(x)\right|=\left|\cos \left(\sqrt{\frac{{\omega }^2}{m{E}_0}\left(\frac{{\rho }_n}{n+1}{x}^{n+1}+{\rho }_{0}x\right)\Phi \left(-\frac{E_m}{E_0}{x}^m;1;\frac{1}{m}\right)}\right)u_0\right|\\ \left|\tau (x)\right|=\left|\sqrt{\frac{{\omega }^2\left(\frac{{\rho }_n}{n+1}{x}^{n+1}+{\rho }_{0}x\right)}{\frac{1}{m{E}_0}\Phi \left(-\frac{E_m}{E_0}{x}^m;1;\frac{1}{m}\right)}}\sin \left(\sqrt{\frac{{\omega }^2}{m{E}_0}\left(\frac{{\rho }_n}{n+1}{x}^{n+1}+{\rho }_{0}x\right)\Phi \left(-\frac{E_m}{E_0}{x}^m;1;\frac{1}{m}\right)}\right)u_0\right|\end{array}$ (3.4)

Аналогично, уравнения удельной энергии (2.27):

$\begin{array}{l}{E}_k=\frac{\pi }{\omega }\left({\rho }_n{x}^n+{\rho }_0\right){\cos }^2\left(\sqrt{\frac{{\omega }^2}{m{E}_0}\left(\frac{{\rho }_n}{n+1}{x}^{n+1}+{\rho }_{0}x\right)\Phi \left(-\frac{E_m}{E_0}{x}^m;1;\frac{1}{m}\right)}\right)u_0^2\\ E_s=\frac{\pi }{\omega \left(E_m{x}^m+E_0\right)}\times \frac{{\omega }^2\left(\frac{{\rho }_n}{n+1}{x}^{n+1}+{\rho }_{0}x\right)}{\frac{1}{m{E}_0}\Phi \left(-\frac{E_m}{E_0}{x}^m;1;\frac{1}{m}\right)}\times \\ \times {\sin }^2\left(\sqrt{\frac{{\omega }^2}{m{E}_0}\left(\frac{{\rho }_n}{n+1}{x}^{n+1}+{\rho }_{0}x\right)\Phi \left(-\frac{E_m}{E_0}{x}^m;1;\frac{1}{m}\right)}\right)u_0^2.\end{array}$ (3.5)

3.2. Линейные биномы. Рассмотрим теперь линейное изменение акустических свойств

$\rho \left(x\right)={\rho }_{1}x+{\rho }_0, E\left(x\right)=E_{1}x+E_0$                    (3.6)

с положительными коэффициентами ${\rho }_0,E_0$ и неотрицательными ${\rho }_1,E_1$. Подстановка линейных членов (3.2) в уравнение (2.18) дает

$a\left(x\right)=\frac{1}{E_0}\text{ln}\left(E_{1}x+E_0\right), b\left(x\right)=-{\omega }^2\left(\frac{{\rho }_1}{2}{x}^2+{\rho }_{0}x\right)$,           (3.7)

откуда получаем

$\left|m\left(x\right)\right|=\left|cos\eta \left(x\right)\right|u_0, \left|\tau \left(x\right)\right|=\left|\sqrt{\frac{{\omega }^2\left(\frac{{\rho }_1}{2}{x}^2+{\rho }_{0}x\right)}{\frac{1}{E_0}\text{ln}\left(E_{1}x+E_0\right)}}sin\eta \left(x\right)\right|u_0$,            (3.8)

где

         $\eta (x)=\sqrt{\frac{{\omega }^2}{E_0}\left(\frac{{\rho }_1}{2}{x}^2+{\rho }_{0}x\right)\ln \left(E_{1}x+E_0\right)}$.              (3.9)

Принимая во внимание уравнения (3.8), где $\eta (x)$ представляет собой аргумент тригонометрической функции, и согласно [41, 42], представим так называемую пространственную частоту

${\Omega }_{\eta }(x)\equiv {\partial }_x\eta (x)=\frac{\sqrt{2}\omega }{4}\frac{2\left(\ln \left(E_{1}x+E_0\right)\right)\left({\rho }_{1}x+{\rho }_0\right)\left(E_{1}x+E_0\right)+x{E}_1\left({\rho }_{1}x+2{\rho }_0\right)}{\sqrt{\left(\ln \left(E_{1}x+E_0\right)\right)x\left({\rho }_{1}x+2{\rho }_0\right)E_0\left(E_{1}x+E_0\right)}}$.   (3.10)

Для предполагаемого случая, выражение для удельной энергии (2.27) представим в виде

$E_k=\frac{\pi }{\omega }\left({\rho }_{1}x+{\rho }_0\right)\left({\text{cos}}^2\eta \right)u_0^2, E_s=\frac{\pi \omega E_0\left(\frac{{\rho }_1}{2}{x}^2+{\rho }_{0}x\right)\left({\text{sin}}^2\eta \right)u_0^2}{\left(E_{1}x+E_0\right)\text{ln}\left(E_{1}x+E_0\right)}$.         (3.11)

Типичные графики вариаций магнитуд и энергии показаны на рис. 1 и 2. Графики на рис. 1 соответствуют возрастающей с расстоянием фазовой скорости $c(x)=\sqrt{E(x)/\rho (x)}$ и следующим безразмерным параметрам материала модели:

${\rho }_0=\mathrm{1,}{ }^{}{\rho }_1=\mathrm{0.005,}{ }^{}{E}_0=\mathrm{1,}{ }^{}{E}_1=\mathrm{0.01.}$                  (3.12)

Графики на рис. 2 показывают, что (i) величина смещения не меняется с расстоянием, что связано с наложенным граничным условием Дирихле на левом конце стержня; (ii) величина напряжения увеличивается (почти линейно в рассматриваемых полулогарифмических координатах), что в основном связано с увеличением модуля упругости; (iii) удельная энергия деформации нелинейно увеличивается с расстоянием, то есть опять же в основном за счет увеличения модуля упругости; и (iv) удельная кинетическая энергия медленно увеличивается в основном за счет увеличения плотности.

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

Графики на рис. 3 соответствуют убывающей с расстоянием фазовой скорости и следующим безразмерным параметрам модельного материала

${\rho }_0=\mathrm{1,}{ }^{}{\rho }_1=\mathrm{0.01,}{ }^{}{E}_0=\mathrm{1,}{ }^{}{E}_1=\mathrm{0.005.}$                  (3.13)

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

4. Заключение. Сравнение величин напряжений и перемещений, показанных для линейных (рис. 2,а, 3,а) биномиальных вариаций соответствующей плотности материала и модуля упругости, показывает, что во всех рассмотренных случаях наблюдаются более выраженные колебания с гораздо меньшей видимой частотой относительно к линейным биномам; это связано с гораздо меньшей пространственной частотой ${\Omega }_{\eta }$, чем ${\Omega }_{\xi }$. В то же время величина смещения остается практически стабильной независимо от биномиала, а величина напряжения демонстрирует относительно небольшой рост с расстоянием. Анализ изменения удельной энергии показывает, что в случаях линейных биномов энергия деформации увеличивается нелинейно, тогда как кинетическая энергия остается практически стабильной (рис. 2,b, 3,b).

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

Полученные результаты могут иметь различные приложения, начиная от микро- и наномеханики при нахождении полей напряжений и смещений и пространственной частоты акустических волн, распространяющихся в микро- и нанотрубках, до геофизических приложений, особенно при количественной оценке неизвестных физических свойств земной коры. Следует также отметить, что, согласно обзору литературы, ни величины полей перемещений, ни пространственное изменение энергии для случая гармонических волн в одномерных ФГ-стержнях не изучались. Однако изменение скорости определяется простой формулой для скорости стержня или, в более общем случае, длинноволновой предельной скоростью, см. [45-47].

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

Работа финансировалась Министерством науки и высшего образования РФ, проект № FSWG-2023-0004 “Система территориальной сейсмической защиты критически важных объектов инфраструктуры на основе гранулированных метаматериалов, обладающих свойствами широкодиапазонных фононных кристаллов”.