Известия РАН. Механика твердого тела, 2022, № 3, стр. 56-69

1. Введение. По современным представлениям, развитым в трудах М.А. Садовского [1] и его последователей, горные породы представляют собой иерархическую систему блоков разного масштабного уровня. Блоки одного уровня разделены между собой прослойками пород с ослабленными механическими свойствами. В [2, 3] отмечено, что размеры блоков изменяются в масштабах от фракций породного массива до гео-блоков земной коры. В экспериментальной работе [4] было показано на двумерной модели блочной среды – кирпичной стене, что для реальной геосреды можно определить размеры характерных блоков породного массива по данным сейсмического каротажа, используя обнаруженное в [5] соотношение, связывающее значение скорости распространения низкочастотной волны, частоты, ограничивающей ее спектр, и продольного размера блоков. Как показано в [2, 3, 6], движение блочной среды может быть представлено, как движение жестких блоков за счет деформации прослоек. В результате динамику блочной среды можно изучать в маятниковом приближении, когда считается, что блоки несжимаемы, а все деформации и смещения происходят за счет сжимаемости прослоек (см., например, [7–9]). В [8, 9] блочная среда моделируется как трехмерная решетка масс, соединенных элементами Фойгта в осевых и диагональных направлениях. В [9] показано качественное соответствие конечно-разностного решения по этой модели задачи Лэмба для блочной среды с результатами полевых экспериментов, проведенных на известняковом карьере. Альтернативный подход основан на математической модели блочной среды с упругими блоками, взаимодействующими через податливые прослойки [5, 10, 11]. Для описания прослоек в [10, 11] предложены различные варианты модели, в которой прослойки между упругими блоками могут быть упругими, вязкоупругими, пластическими, пористыми.

В настоящей работе модифицируется 3D модель, предложенная в [8]. В новой модели внутреннее трение в прослойках между блоками моделируется элементами Максвелла и Фойгта. Зенер [12] и Био [13] были одними из первых, кто включили элементы Максвелла и Фойгта в модели для описания неупругого поведения материалов. Модель с двумя элементами Максвелла и одним элементом Фойгта впервые была предложена Био [13], кроме того она была представлена в работе Фанг [14]. Многие исследователи затем использовали различные комбинации элементов Максвелла и Фойгта для учета неупругих потерь. При этом важной задачей в течение нескольких десятилетий в каждой модели было ограничение количества механизмов релаксации для получения удовлетворительного решения, т.е. для получения почти постоянного коэффициента добротности материала Q. В [15] исследована модель с внутренним трением для описания распространения волн в однородных неупругих средах. Эта модель включает два элемента Максвелла и один элемент Фойгта. Как показано в [15], эта модель дает возможность выбирать параметры вязкости и жесткости этих элементов, так что коэффициент добротности материала Q отличается от предписанного постоянного значения ${{Q}_{0}}$ не больше, чем на 5% на интервале частот от 3% до 100% от максимальной частоты, представляющей интерес. Важным свойством данной модели является небольшое число дополнительных переменных при максимальном охвате интересующего диапазона частот. Использование минимального количества элементов Максвелла ограничивает память и вычислительные ресурсы, необходимые в компьютерных программах для численного моделирования трехмерных задач.

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

2. Одномерная модель блочной среды с учетом внутреннего трения. Продемонстрируем эту модель сначала на примере дискретно-периодической одномерной цепочки масс, соединенных вязкоупругими прослойками. Реологическая модель прослоек состоит из двух элементов Максвелла и одного элемента Фойгта (рис. 1). Все три элемента соединены параллельно. Каждый элемент Максвелла состоит из пружины и демпфера, которые соединены последовательно. Элемент Фойгта состоит также из пружины и демпфера, но соединены они параллельно. Здесь k, k₁, k₂ – жесткости пружин в элементе Фойгта и в двух элементах Максвелла, соответственно, $\lambda $, ${{\lambda }_{1}}$, ${{\lambda }_{2}}$ – вязкости демпферов в элементе Фойгта и в двух элементах Максвелла.

Уравнения одномерного движения блоков с данной реологической моделью прослоек имеют следующий вид:

Здесь u_j – перемещения жестких блоков, $P(t)$ – действующая нагрузка, приложенная к блоку $j = 0$; K – суммарная жесткость пружин. В (2.1) введены две дополнительные переменные ${{\varphi }_{j}}$ и ${{\psi }_{j}}$, которые зависят от функции смещения u_j.

Применим преобразование Фурье по времени к уравнениям движения (2.1). В частотной области ω сила, действующая на блок со стороны прослойки, может быть выражена формулой:

где есть нормированная функция импеданса [16].

Внутреннее затухание в среде обычно определяется коэффициентом добротности Q(ω) или его обратной величиной Q⁻¹(ω), которая описывается формулой:

Для почв и породных материалов принято предполагать, что целевой коэффициент добротности материала остается постоянным в широком диапазоне интересующих частот, т.е. ${{Q}_{0}}(\omega ) = {{Q}_{0}}$, см. [17].

Чтобы связать параметры ${{\alpha }_{1}}$, ${{\alpha }_{2}}$, ${{\gamma }_{1}}$, ${{\gamma }_{2}}$ и $\beta $ реологической модели с Q⁻¹(ω) с помощью простых, не зависящих от частоты приближений, введем новые безразмерные параметры ${{\hat {\alpha }}_{1}}$, ${{\hat {\alpha }}_{2}}$, ${{\hat {\gamma }}_{1}}$, ${{\hat {\gamma }}_{2}}$, $\hat {\beta }$ и частоту $\hat {\omega }$ по формулам:

где ${{\omega }_{{\max }}}$ – максимальная угловая частота, представляющая интерес для моделирования.

В выражении (2.2) коэффициента ${{Q}^{{ - 1}}}(\omega )$ перейдем к нормированным параметрам (2.3). Тогда формула (2.2) примет вид:

Из этой формулы можно увидеть, что коэффициенты ${{\hat {\alpha }}_{1}}$, ${{\hat {\alpha }}_{2}}$, ${{\hat {\gamma }}_{1}}$, ${{\hat {\gamma }}_{2}}$, $\hat {\beta }$ не зависят от ${{\omega }_{{\max }}}$. Отсюда следует, что для заданной постоянной $Q_{0}^{{ - 1}}$, они должны вычисляться только один раз. При больших значениях целевого коэффициента добротности, то есть $Q_{0}^{{ - 1}} \to 0$, параметры становятся существенно не зависимыми от Q₀.

Задача определения параметров ${{\hat {\alpha }}_{1}}$, ${{\hat {\alpha }}_{2}}$, ${{\hat {\gamma }}_{1}}$, ${{\hat {\gamma }}_{2}}$, $\hat {\beta }$ решается с допуском, так чтобы фактический коэффициент добротности оставался бы достаточно близким к целевому значению. Этот допуск, который принимается как 5% от целевого значения, выполняется в пределах от 4% до 100% от ω_max. Т.е. для заданного коэффициента добротности ${{Q}_{0}}$ параметры ${{\hat {\alpha }}_{1}}$, ${{\hat {\alpha }}_{2}}$, ${{\hat {\gamma }}_{1}}$, ${{\hat {\gamma }}_{2}}$, $\hat {\beta }$ подбираются прямой проверкой на компьютере так, чтобы условие

выполнялось для каждого $\hat {\omega }$, удовлетворяющего неравенству $0.04 \leqslant \hat {\omega } \leqslant 1$. Найденные значения параметров, соответствующие пяти значениям целевого коэффициента добротности $Q_{0}^{{ - 1}}$, представлены в табл. 1.

Для этих значений параметров фактические коэффициенты добротности материала, нормированные относительно $Q_{0}^{{ - 1}}$, показаны как функции частоты на рис. 2. Как видно на рис. 2, удается подобрать параметры ${{\hat {\alpha }}_{1}}$, ${{\hat {\alpha }}_{2}}$, ${{\hat {\gamma }}_{1}}$, ${{\hat {\gamma }}_{2}}$, $\hat {\beta }$ так, что допуск, который принимается как 5% от целевого значения, выполняется в пределах от 4% до 100% от ${{\omega }_{{\max }}}$.

Другим важным аспектом внутренних моделей трения является дисперсия, наблюдаемая во временном отклике системы. Для дискретного уравнения (2.1) дисперсионная зависимость волнового числа q от частоты ω имеет вид:

Здесь l – длина пружин, q – волновое число. Пользуясь формулой (2.4) и соотношением ${{c}_{{{\text{ph}}}}}(\omega ) = {\omega \mathord{\left/ {\vphantom {\omega q}} \right. \kern-0em} q}$, где ${{c}_{{{\text{ph}}}}}$ – фазовая скорость, можем определить зависимость фазовой скорости от частоты для одномерной модели блочной среды (2.1). Если положить $\omega \to 0$, то из (2.4) получим: ${{c}_{{{\text{ph}}}}}(0) = l\sqrt {{{K{\kern 1pt} \left( {1 - {{\alpha }_{1}} - {{\alpha }_{2}}} \right)} \mathord{\left/ {\vphantom {{K{\kern 1pt} \left( {1 - {{\alpha }_{1}} - {{\alpha }_{2}}} \right)} M}} \right. \kern-0em} M}} $. Как следует из этой формулы и таблицы 1, фазовая скорость бесконечно длинных волн падает с ростом $Q_{0}^{{ - 1}}$.

На рис. 3 представлены зависимости ${{c}_{{{\text{ph}}}}}(\omega )$ при различных значениях коэффициента добротности материала $Q_{0}^{{ - 1}}$. Здесь и далее полагаем: $l = 1$, $M = 1$, $K = 1$. В прямоугольнике представлены эти же кривые в увеличенном масштабе в низкочастотной зоне. Заметим, что в дискретном случае имеем: ${{\omega }_{{\max }}} = 2{{\omega }_{0}}$, где ${{\omega }_{0}} = \sqrt {{K \mathord{\left/ {\vphantom {K M}} \right. \kern-0em} M}} $. Видно, что наибольшее влияние коэффициента добротности на дисперсию наблюдается в низкочастотной части спектра. Кроме того, видно, что с ростом $Q_{0}^{{ - 1}}$ максимальная скорость распространения возмущений падает. В высокочастотной части спектра фазовая скорость практически не зависит от коэффициента добротности.

На рис. 4 представлены результаты численных расчетов нестационарного волнового процесса в цепочке масс при действии нагрузки $P(t) = H(t)$, где $H(t)$ – ступенчатая функция Хевисайда. Расчеты проведены по одномерной модели с внутренним трением (2.1) для блока с номером $j = 10$. Сплошные кривые соответствуют $Q_{0}^{{ - 1}} = 0$, штрих-пунктирные кривые – $Q_{0}^{{ - 1}} = 0.2$. На рис. 4, a показаны деформации прослойки $\Delta {{u}_{j}} = {{u}_{{j + 1}}}$ – u_j, на рис. 4, b – скорости ${{\dot {u}}_{j}}$ и на рис. 4, c – ускорения ${{\ddot {u}}_{j}}$ в зависимости от времени.

Горизонтальные линии соответствуют статическим значениям: на рис. 4, a – $\Delta {{u}_{{{\text{st}}}}} = - {1 \mathord{\left/ {\vphantom {1 {\left( {1 - {{\alpha }_{1}} - {{\alpha }_{2}}} \right)}}} \right. \kern-0em} {\left( {1 - {{\alpha }_{1}} - {{\alpha }_{2}}} \right)}}$, на рис. 4, b – ${{\dot {u}}_{{{\text{st}}}}} = {1 \mathord{\left/ {\vphantom {1 {\sqrt {1 - {{\alpha }_{1}} - {{\alpha }_{2}}} }}} \right. \kern-0em} {\sqrt {1 - {{\alpha }_{1}} - {{\alpha }_{2}}} }}$. К этим значениям стремятся амплитуды деформации прослоек и скоростей блоков при больших значениях времени. Видно, что с ростом $Q_{0}^{{ - 1}}$ увеличивается затухание максимальных амплитуд деформаций прослоек, скоростей и ускорений блоков на фронте низкочастотной волны и падает амплитуда высокочастотных осцилляций $\Delta {{u}_{j}}$ и ${{\dot {u}}_{j}}$ относительно их статических значений за фронтом волны.

3. Трехмерная модель блочной среды с учетом внутреннего трения. Блочная среда моделируется однородной трехмерной решеткой, состоящей из точечных масс, соединенных пружинами и демпферами в направлениях осей x, y, z и в диагональных направлениях плоскостей xy, xz, yz, как показано на рис. 5, a. Здесь использованы обозначения: u, ${v}$ – горизонтальные перемещения в направлениях x, y; w – вертикальные перемещения в направлении z; n, m, k – номера блоков в направлениях x, y, z. На поверхности блочного полупространства в начале координат приложена вертикальная сосредоточенная нагрузка P (рис. 5, b).

В качестве реологической модели прослоек используется модель прослоек с двумя элементами Максвелла и одним элементом Фойгта, которая только что была продемонстрирована на примере одномерной блочной среды. Далее будем полагать, что жесткости пружин и вязкости демпферов в осевых и диагональных направлениях совпадают. Кроме того, будем полагать, что параметры M, k, ${{k}_{1}}$, ${{k}_{2}}$, $\lambda $, ${{\lambda }_{1}}$, ${{\lambda }_{2}}$ имеют одни и те же значения во всех точках блочной среды.

С учетом внутреннего трения и обозначений

уравнения движения блока с координатами n, m, k, находящегося внутри полупространства $(k < 0)$, имеют вид:

С учетом обозначений

уравнения движения блока с координатами n, m, 0, находящегося на свободной поверхности полупространства, имеют вид:

Здесь ${{\delta }_{{n0}}}$ – символ Кронекера, $P(t)$ – амплитуда, действующей поверхностной нагрузки. Начальные условия нулевые.

Как показано в [8], уравнениям (3.1) – (3.6), описывающим движение блочного полупространства с упругими прослойками (${{\alpha }_{1}} = 0$, ${{\alpha }_{2}} = 0$, ${{\gamma }_{1}} = 0$, ${{\gamma }_{2}} = 0$, $\beta = 0$), соответствуют скорости продольных (${{c}_{{\text{p}}}}$), сдвиговых (${{c}_{{\text{s}}}}$) и рэлеевских (${{c}_{{\text{R}}}}$) бесконечно длинных волн:

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

Для этих задач исследуется реакция блочной среды на два вида действующей нагрузки: (а) ступенчатое воздействие, т.е. $P(t) = {{P}_{0}}H(t)$; (б) импульс Гаусса, т.е. P(t) = = ${{P}_{0}}{\text{exp}}[ - {{{{{(t - 4\sigma )}}^{2}}} \mathord{\left/ {\vphantom {{{{{(t - 4\sigma )}}^{2}}} {(2{{\sigma }^{2}})}}} \right. \kern-0em} {(2{{\sigma }^{2}})}}]$.

4. Разностная схема. Уравнения (3.1)–(3.6) с нулевыми начальными данными решались методом конечных разностей по явной схеме. Для вторых производных по времени использовалась центрально-разностная аппроксимация второго порядка точности, для первых производных по времени применялась разность “назад” первого порядка точности:

Здесь $\tau $ – шаг разностной сетки по времени, $u_{{n,m,k}}^{s}$ – значение перемещения ${{u}_{{n,m,k}}}(t)$ в момент времени $t = s\tau $, s – номер слоя по времени в конечно-разностной схеме. Для дополнительных функций использовались аппроксимации следующего вида:

5. Результаты численных расчетов задачи 1 с учетом внутреннего трения. В этой секции представлены результаты численных расчетов возмущений в задаче Лэмба с учетом внутреннего трения, проведенных при ${{P}_{0}} = 1$ и $\tau = \pi {\text{/}}20$. Масса блоков, длина пружин и суммарная жесткость прослоек приняты за единицы: $M = 1$, $l = {\text{1}}$, $K = 1$.

На рис. 6 – 9 приведены результаты расчетов вертикальных $\dot {w} = {{\dot {w}}_{{n,m,k}}}$ и радиальных ${{\dot {u}}_{r}} = {{\left( {n\,{{{\dot {u}}}_{{n,m,k}}} + m\,{{{{\dot {v}}}}_{{n,m,k}}}} \right)} \mathord{\left/ {\vphantom {{\left( {n\,{{{\dot {u}}}_{{n,m,k}}} + m\,{{{{\dot {v}}}}_{{n,m,k}}}} \right)} {{{{({{n}^{2}} + {{m}^{2}})}}^{{1/2}}}}}} \right. \kern-0em} {{{{({{n}^{2}} + {{m}^{2}})}}^{{1/2}}}}}$ скоростей блоков в цилиндрической системе координат $r,\;\theta ,\;z$ на поверхности блочного полупространства (k = 0) на оси $m = 0$.

На рис. 6 приведены графики зависимости от времени радиальной ${{\dot {u}}_{r}}$ и вертикальной $\dot {w}$ скоростей блоков в точке (60, 0, 0), рассчитанные при действии ступенчатой нагрузки при различных значениях коэффициента добротности материала $Q_{0}^{{ - 1}}$. Вертикальные линии соответствуют моментам времени прихода продольных волн ${{t}_{{\text{p}}}} = {n \mathord{\left/ {\vphantom {n {{{c}_{{\text{p}}}}}}} \right. \kern-0em} {{{c}_{{\text{p}}}}}}$ и рэлеевских волн ${{t}_{{\text{R}}}} = {n \mathord{\left/ {\vphantom {n {{{c}_{{\text{R}}}}}}} \right. \kern-0em} {{{c}_{{\text{R}}}}}}$, где ${{c}_{{\text{p}}}}$ и ${{c}_{{\text{R}}}}$ определяются по формулам (3.7). Видно, что с ростом $Q_{0}^{{ - 1}}$ увеличивается затухание амплитуды скоростей блоков для всего спектра частот.

На рис. 7 приведены результаты расчетов максимальных значений абсолютных величин амплитуд радиальной и вертикальной скоростей блоков при действии ступенчатой по времени нагрузки. На рис. 7, a, b показаны зависимости этих величин от координаты n для различных значений коэффициента добротности $Q_{0}^{{ - 1}}$, на рис. 7, c – зависимости от коэффициента добротности $Q_{0}^{{ - 1}}$ в точке (90, 0, 0). Видно, что при увеличении коэффициента $Q_{0}^{{ - 1}}$ максимальные значения абсолютных величин амплитуд ${{\dot {u}}_{r}}$ и $\dot {w}$ падают быстрее с ростом n и их зависимость от коэффициента $Q_{0}^{{ - 1}}$ близка к линейной.

На рис. 8 приведены графики зависимости от времени радиальной и вертикальной скоростей блоков в точке (60, 0, 0), рассчитанные при различных значениях коэффициента добротности $Q_{0}^{{ - 1}}$ при действии импульса Гаусса ($\sigma = 5$). Вертикальные линии соответствуют моментам прихода продольных и рэлеевских волн. Видно, что с ростом $Q_{0}^{{ - 1}}$ уменьшаются амплитуды скоростей блоков и увеличивается время наступления максимальных амплитуд скоростей. Последнее объясняется уменьшением скорости низкочастотных волн, что можно видеть на рис. 2.

На рис. 9 приведены результаты расчетов максимальных значений абсолютных величин амплитуд радиальной и вертикальной скоростей блоков на поверхности блочного полупространства при действии импульса Гаусса ($\sigma = 5$). На рис. 9, a, b показаны зависимости этих величин от координаты n для различных значений коэффициента добротности $Q_{0}^{{ - 1}}$, на рис. 9, c – зависимости от коэффициента добротности $Q_{0}^{{ - 1}}$ в точке (100, 0, 0).

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

Анализ этих приближенных функций и графиков на рис. 9, a, b показывает, что с ростом коэффициента $Q_{0}^{{ - 1}}$ скорость затухания максимальных по модулю значений ${{\dot {u}}_{r}}$, $\dot {w}$ с удалением от места воздействия растет. На рис. 9, c видно, что зависимость максимальных по модулю амплитуд скоростей блоков от коэффициента $Q_{0}^{{ - 1}}$ близка к параболической.

6. Результаты численных расчетов задач 1 и 2 с учетом внутреннего трения. Ниже, на рис. 10 представлены результаты численных расчетов двух задач с учетом внутреннего трения при воздействии вертикального сосредоточенного импульса Гаусса, во-первых, на поверхность блочного полупространства (задача 1 – штриховые кривые) и, во-вторых, на поверхность блочного слоя, лежащего на блочном полупространстве, имеющем другие свойства (задача 2 – сплошные кривые).

На рис. 10 показаны зависимости от времени радиальных и вертикальных скоростей блоков на поверхности блочного полупространства, находящихся на различном расстоянии от места воздействия (n = 20, 40, 60, 80). Коэффициент добротности $Q_{0}^{{ - 1}}$ = 0.1 для блочного полупространства (задача 1) и для слоя, лежащего на полупространстве (задача 2). Коэффициент добротности $Q_{0}^{{ - 1}} = 0.05$ для блочного полупространства, лежащего под слоем (задача 2). Кроме того, скорость распространения волн в блочном полупространстве и в блочном слое, лежащем на полупространстве, в два раза меньше, чем скорость распространения волн в полупространстве, лежащем под слоем.

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

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

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

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

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

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

Работа выполнена в рамках проекта НИР, номер государственной регистрации 121062200075-4.