Акустический журнал, 2020, T. 66, № 3, стр. 235-241

ВВЕДЕНИЕ

Задачи о колебаниях упругих полуограниченных тел [1–11] достаточно давно привлекают внимание многих исследователей, поскольку анализ динамических процессов весьма важен, например, при оценке динамической прочности элементов тяжелонагруженных транспортных средств, машин и механизмов, в геофизике и горной механике. Отметим, что для изотропных однородных сред динамические процессы в таких структурах хорошо изучены [1, 2], а исследование волновых процессов для сред с более сложной структурой, например, для слоистых [3, 4], функционально-градиентных [5, 6], пористых [7–9] является весьма актуальным и в настоящее время. Одним из важных классов динамических задач теории упругости являются задачи о колебаниях упругих тел с локальными неоднородностями [10–13], такими как полости, включения, трещины. Изучение рассеяния упругих волн на неоднородностях разных типов и форм позволяет устанавливать факт наличия дефекта [14–16], что дает возможность определять тип и положение повреждения в исследуемой структуре. Чаще всего решение таких задач базируется на сведении их к системам граничных интегральных уравнений [1, 2, 11, 12] при помощи фундаментальных решений или матриц-функций Грина, что позволяет снизить размерность исходной задачи на единицу. Такой подход к решению реализован в настоящей работе. При решении динамических задач теории упругости также применяется метод блочных элементов [3] и гибридные численно-аналитические схемы [10, 13]. На основе решения задач о колебаниях сред с дефектами может быть осуществлен расчет амплитудных значений волновых полей, используемых для идентификации дефекта. В настоящей работе реализован подход, основанный на расчете волнового поля в волноводе с полостью на основе метода граничных элементов [17, 18], асимптотического подхода и приближения Борна, часто используемого в акустике [19, 20], и проведено сравнение результатов.

ПОСТАНОВКА И РЕШЕНИЕ ЗАДАЧИ

Рассмотрим установившийся режим колебаний с частотой $\omega $ ортотропного упругого слоя толщины $h$ с цилиндрической полостью, не выходящей на его границы, направляющая которой есть гладкая замкнутая кривая ${{l}_{0}}$ с образующей, параллельной оси $O{{x}_{2}}$. Нижняя грань слоя жестко защемлена и совпадает с осью $О{{х}_{1}}$, ось $О{{х}_{3}}$ направлена перпендикулярно вверх. Оси упругой симметрии материала совпадают с осями системы координат. Колебания в слое вызваны нагрузкой ${\mathbf{p}}({{x}_{1}},t) = \operatorname{Re} ({\mathbf{p}}({{x}_{1}}){{e}^{{ - i\omega t}}}),$ приложенной к верхней части границы слоя.

Краевая задача после отделения временного множителя имеет вид:

где $\rho $ – плотность среды, ${{C}_{{ijkl}}}$ – компоненты тензора упругих постоянных материала, удовлетворяющие соотношениям симметрии и положительной определенности, ${{n}_{j}}$ – компоненты единичного вектора нормали к кривой ${{l}_{0}}$, внешнего по отношению к области, занятой упругой средой. Замыкает постановку задачи условие излучения волн на бесконечности, при формулировке которого использован принцип предельного поглощения [2].

Исходная задача (1)–(2) в зависимости от типа нагрузки распадается на две: задача об антиплоских колебаниях ортотропного слоя с цилиндрической полостью (задача 1) и плоская задача о колебаниях слоя с полостью (задача 2). В задаче 1 ненулевой является компонента ${{u}_{2}}({{x}_{1}},{{x}_{3}})$ и в краевой задаче (1)–(2) полагаем $i = 2,$ $j = 1,3$. Для задачи 2 отличными от нуля являются компоненты ${{u}_{i}}({{x}_{1}},{{x}_{3}}),$ $i = 1,3$ и в краевой задаче (1)–(2) полагаем $i,j = 1,3.$

Основным способом исследования задач 1, 2 является предварительное сведение их к интегральным уравнениям с нерегулярными ядрами на основе идей теории потенциала. Решения краевых задач 1, 2 построены на основе матриц-функций Грина для слоя $U_{i}^{{(m)}}(x,\xi ),\,\,\,\,i = 1,2,3$, удовлетворяющих однородным граничным условиям на границах, и обобщенной теоремы взаимности [2], где индекс $m$ указывает, что сосредоточенная сила приложена в направлении оси $O{{x}_{m}},\,\,m = 1,2,3$. Функция Грина для ортотропного слоя в случае задачи 1 приведена в работе [21], а для задачи 2 имеет вид:

Здесь $\Delta _{i}^{{(m)}},\,\,i,\,\,m = 1,3$ – определители, получающиеся путем замены $\,\,k$-го столбца ($k = 1,2$) в (4) столбцом ${{g}^{{(m)}}} = {{\left( {{{\delta }_{{1m}}}{{\delta }_{{3m}}}} \right)}^{T}},$ ${{\alpha }_{i}},\,\,i = 1,3$ – параметры преобразования Фурье, $\Gamma $ – поверхность, всюду совпадающая с плоскостью ${{R}^{2}}$, за исключением множества нулей $\Delta $, которые она огибает в соответствии с принципом предельного поглощения. Второе слагаемое в выражении (3) представляет собой некоторую регулярную при $ - h + {{\delta }_{1}} < {{x}_{3}},$ ${{\xi }_{3}} < h - {{\delta }_{1}},$ ${{\delta }_{1}} > 0$ добавку, обеспечивающую выполнение однородных граничных условий. Надо отметить, что представление функций Грина для ортотропного слоя в виде однократного интеграла по контуру $\sigma $ в комплексной плоскости [21, 22] удобно использовать при построении решения задачи (1)–(2) путем сведения ее к интегральным уравнениям с нерегулярными ядрами (6), тогда как при решении задачи в рамках асимптотического подхода использовался вид (3).

Поле перемещений в слое под действием поверхностной нагрузки с носителем на отрезке $[a,b]$ для задач 1, 2 имеет вид:

где $\sigma _{{ij}}^{{(m)}}({\mathbf{x}},\xi )$ – компоненты тензора напряжений, выраженные через функции Грина для слоя $U_{i}^{{(m)}}({\mathbf{x}},\xi ),$ $i = 1,2,3.$

Поля перемещений в слое (5) представимы в виде суммы двух слагаемых, первое из которых $u_{m}^{*}(\xi )$ – поле смещений в среде без дефекта под действием заданной нагрузки (эталонное поле смещений), второе слагаемое обусловлено наличием полости в слое. Таким образом, на основании соотношения (5) может быть рассчитано поле смещений всюду в слое, если найдено поле смещений на границе полости. В настоящей работе поля смещений на контуре полости рассчитаны на основании двух подходов: метод интегральных уравнений и основанный на нем метод граничных элементов (МГЭ), который может быть использован для полостей любой конфигурации, и асимптотический метод для круговых полостей малого относительного радиуса.

Далее сформулированы системы граничных интегральных уравнений (ГИУ), которые за счет специального выбора матриц-функций Грина содержат операторы лишь по границе полости. Системы ГИУ для задачи 1 $(j = 1,3;\,\,i{\text{,}}m = 2)$ и задачи 2 ($i,j,m = 1,3$) имеют вид:

где интеграл по ${{l}_{0}}$ понимается в смысле главного значения по Коши, а в случае нерегулярной кривой ${{l}_{0}}$ коэффициент при ${{u}_{m}}({\mathbf{y}})$ отличен от ½ и зависит от угла между касательными в особой точке и упругих постоянных материала [22].

В настоящей работе численная реализация систем ГИУ (6) с нерегулярным ядром на основе МГЭ [17, 18] выполнена следующим образом. Гладкая граница аппроксимировалась N-угольником $l = \bigcup\nolimits_{q = 1}^N {{{l}_{q}}} $, где функция ${{u}_{m}}({{x}_{1}},{{x}_{3}})$ постоянна в пределах каждого из элементов ${{l}_{q}}$, а в качестве узлов выбраны середины соответствующих отрезков. Число ГЭ $N$ выбиралось таким образом, чтобы было не менее 5–6 элементов на длину волны. Таким образом, системы ГИУ (6) сводились к линейным алгебраическим системам относительно ${{u}_{{mq}}} = {{u}_{m}}({{y}_{{1q}}},{{y}_{{3q}}}).$ Полученные алгебраические системы относительно узловых смещений на контуре хорошо обусловлены, имеют явное диагональное преобладание. Отметим, что МГЭ позволяет рассчитывать поле перемещений на поверхности слоя в случае произвольной формы и размера дефекта. Процедура решения ГИУ (6) достаточно сложна и требует значительной вычислительной работы, однако, если сечение полости представляет собой окружность малого характерного размера, то схему расчета волновых полей можно значительно упростить на основе асимптотического анализа.

Асимптотический подход к решению многопараметрической задачи (1)–(2) о колебаниях ортотропного слоя с цилиндрической полостью, поперечное сечение которой представляет собой окружность радиуса $r$ с центром в точке ${{{\mathbf{x}}}_{0}} = ({{x}_{{10}}},{{x}_{{30}}}),$ был реализован в области изменения безразмерных параметров:

где ${{\varepsilon }_{1}} = {r \mathord{\left/ {\vphantom {r h}} \right. \kern-0em} h},$ ${{\varepsilon }_{2}} = \omega r\sqrt {{\rho \mathord{\left/ {\vphantom {\rho C}} \right. \kern-0em} C}} ,$ где $C = {{C}_{{44}}}$ для задачи 1, $C = {{C}_{{33}}}$ для задачи 2, а $\varepsilon *$ соответствует частоте запирания, когда в слое нет бегущих волн.

Поля перемещений, всюду в области, занятой упругой средой, отличные от соотношений (5) и полученные на основании обобщенной теоремы взаимности, имеют вид:

где ${{S}_{0}}$ – плоская область, ограниченная гладкой замкнутой кривой ${{l}_{0}}$. На основании представления (8) при предельном переходе $\xi \to {\mathbf{y}} \in {{l}_{0}}$ получена следующая система ГИУ:

Надо отметить, что система ГИУ (9) имеет регулярные подынтегральные функции в отличие от уравнений (6), а криволинейный интеграл по контуру ${{l}_{0}}$ существует как несобственный, что позволяет при дискретизации избежать вычисления сингулярных интегралов. Функцию ${{K}_{{im}}}({\mathbf{x}},{\mathbf{y}})$ можно представить в виде суммы двух слагаемых, первое из которых соответствует статическому случаю (${{\varepsilon }_{2}} = 0$), а второе – некоторая добавка:

При этом имеют место следующие оценки:

Выше введены следующие обозначения для задачи 1:

а для задачи 2:

Поля смещений на контуре полости ${{l}_{0}}\,$ представимы в виде разложения по малому параметру $\,{{\varepsilon }_{1}}$:

Коэффициенты в разложении (10) вычислены аналитически и для задачи 1 приведены в работах [11, 12], а для задачи 2 имеют вид:

Интеграл $I({{C}_{1}},{{C}_{2}},{{C}_{3}}) = \frac{{{{C}_{1}}{{\Psi }_{1}} + {{C}_{2}}{{\Psi }_{2}} + {{C}_{3}}{{\Psi }_{3}}}}{\Lambda }$ вычислен аналитически по теории вычетов, где

Таким образом, для задач 1, 2 получено явное асимптотическое представление полей смещений на контуре полости, минуя процедуру дискретизации ГИУ (6) на основе МГЭ. Поля перемещений на поверхности ${{u}_{m}}({{\xi }_{1}},h),$ $m = 1,2,3$, построенные в рамках асимптотического подхода, в дальней от дефекта зоне (${{x}_{1}} < {{\xi }_{1}}$), имеют вид:

где $M$ – число распространяющихся волн в слое, $\alpha _{{1p}}^{*}$ – корни дисперсионного уравнения. Амплитуды бегущих волн на поверхности слоя для задачи 1 имеют вид [12], а для задачи 2 имеют вид:

Выражения $\psi (\alpha _{{1p}}^{*}),$ $\Phi _{k}^{{(m)}}({{x}_{{30}}}),$ $\Omega _{k}^{{(m)}}({{x}_{{30}}}),$ $\Lambda _{k}^{{(m)}}({{x}_{{30}}})$ не приводятся в виду их громоздкости, а ${{u}_{{s,k}}}({{{\mathbf{x}}}_{0}}),$ $s,k = 1,3$ заданы формулами (11).

Надо отметить, что в рамках приближения Борна, широко используемого в акустике [19, 20], перемещения на контуре заменяются эталонными полями смещений в точке ${{{\mathbf{x}}}_{0}}$, поэтому амплитуды полей перемещений на поверхности слоя в этом приближении не учитывают изменяемость полей смещений на контуре. В рамках такого приближения однократного рассеяния поля смещений на контуре полости определяются соотношением (10) при удержании лишь первого члена в разложении полей. В свою очередь, в соотношении (13) для амплитуд отсутствуют слагаемые, содержащие производные от эталонных полей $u_{m}^{*}({{{\mathbf{x}}}_{0}}),\,\,m = 1,3,$ а их выражение имеет вид:

ЧИСЛЕННЫЙ ЭКСПЕРИМЕНТ

На основании соотношений (5), (12)–(14) построены поля перемещений ${{\upsilon }_{m}}({{\xi }_{1}},h)$ = ${{u}_{m}}({{\xi }_{1}},h) - u_{m}^{*}({{\xi }_{1}},h),$ $m = 1,2,3$ на поверхности слоя ${{x}_{3}} = h,$ обусловленные лишь наличием полости. Проведен ряд численных экспериментов по расчету волновых полей на поверхности слоя для случая цилиндрической полости кругового сечения малого относительного размера на основании трех походов: МГЭ, асимптотического подхода и приближения Борна.

В расчетах приняты следующие значения упругих постоянных для диоксида гафния HfО₂ ${{С}_{{11}}} = 502,$ ${{С}_{{13}}} = 159,$ ${{С}_{{33}}} = 597,$ ${{C}_{{44}}} = 78,$ ${{С}_{{55}}} = 90,$ ${{C}_{{66}}} = 111$ (ГПа), $\rho = 13\,310$ кг/м³ [22]. Колебания в слое вызваны сосредоточенной нагрузкой ${\mathbf{p}} = \{ 0,{{ - {{p}_{0}}} \mathord{\left/ {\vphantom {{ - {{p}_{0}}} {{{C}_{{44}}},}}} \right. \kern-0em} {{{C}_{{44}}},}}0\} $ и ${\mathbf{p}} = \{ 0,0,{{ - {{p}_{0}}} \mathord{\left/ {\vphantom {{ - {{p}_{0}}} {{{C}_{{33}}}}}} \right. \kern-0em} {{{C}_{{33}}}}}\} $ для задач 1 и 2 соответственно, которая приложена в точке $(0,h)$. Направляющая цилиндрической полости ${{l}_{0}}$ – окружность с центром в точке $(0,\,\,{h \mathord{\left/ {\vphantom {h 2}} \right. \kern-0em} 2}).$ Число граничных элементов $N = 12$. На рис. 1–4 приведены графики полей перемещений на поверхности слоя ${{\upsilon }_{m}}({{\xi }_{1}},h),$ $m = 1,2,3$, рассчитанных тремя способами для задач 1 и 2, при этом сплошная линия соответствует расчетам по МГЭ, графики полей, полученных на основе асимптотического подхода, изображены сплошной линией с маркером “*” и штриховой линией отмечены поля, полученные в рамках приближения Борна. На рис. 1–2 представлены графики, соответствующие случаю одной распространяющейся моды в слое для значений параметров ${{\varepsilon }_{1}} = 0.24,$ ${{\varepsilon }_{2}} = 0.98$ (задача 1) и ${{\varepsilon }_{1}} = 0.15,$ ${{\varepsilon }_{2}} = 0.3$ (задача 2).

На рис. 3–4 представлены графики полей перемещений на поверхности слоя ${{\upsilon }_{m}}({{\xi }_{1}},h),$ $m = 1,2,3$ для двух распространяющихся мод при следующих значениях параметров ${{\varepsilon }_{1}} = 0.15,$ ${{\varepsilon }_{2}} = 0.98$ (задача 1) и ${{\varepsilon }_{1}} = 0.1,$ ${{\varepsilon }_{2}} = 0.4$ (задача 2).

ВЫВОДЫ

На основе численных экспериментов определена область корректной работы асимптотического подхода для задач 1 и 2 для среднезаглубленных полостей малого поперечного размера и приближения Борна при расчете волновых полей в дальней зоне. Эта область описывается следующими неравенствами: $\left\{ \begin{gathered} {{\varepsilon }_{1}} < 0.24 \hfill \\ {{\varepsilon }_{2}} < 1 \hfill \\ \end{gathered} \right.$ для задачи 1 и $\left\{ \begin{gathered} {{\varepsilon }_{1}} < 0.15 \hfill \\ {{\varepsilon }_{2}} < 1 \hfill \\ \end{gathered} \right.$ для задачи 2 (относительное расхождение результатов не превышало 10%). Как видно из графиков на рис. 1, на частотах, соответствующих одной распространяющейся моде в слое, волновое поле в дальней от дефекта зоне, рассчитанное на основе приближения Борна, практически совпадает с полем, рассчитанным на основе МГЭ. На тех частотах, для которых имеются две распространяющиеся моды, приближение Борна значительно отличается от расчетов по МГЭ. В то же время поле на поверхности, рассчитанное на основе асимптотического подхода, отличается от поля, полученного на основе МГЭ, не более, чем на 4%. Отметим, что для задачи 2 для любых частот, соответствующих распространяющим модам в слое, наблюдается качественное и количественное отличие полей перемещений, построенных на основании приближения Борна от полей, рассчитанных на основе МГЭ и асимптотического подхода.

Работа выполнена в рамках федеральной целевой программы “Исследования и разработки по приоритетным направлениям развития научно-технологического комплекса России на 2014–2020 годы” при финансовой поддержке государства в лице Министерства науки и высшего образования России (идентификатор проекта RFMEFI60718X0203).