CC BY
22
УДК 539.375
Модель динамического нагружения полости в трехмерной
блочной среде
Н.И. Александрова
Институт горного дела им. Н.А. Чинакала СО РАН, Новосибирск, 630091, Россия
Изучается процесс распространения волн при динамическом нагружении поверхности полости в блочной среде. Предложена трехмерная математическая модель нестационарного вязкоупругого деформирования блочной среды. Эта модель основана на представлении, что динамическое поведение блочной среды может быть приближенно описано как движение жестких блоков за счет податливости прослоек между ними и что деформацию прослоек можно приближенно описать моделью Кельвина-Фойгта. В результате блочная среда моделируется однородной трехмерной решеткой, состоящей из точечных масс, соединенных пружинами и демпферами в направлениях осей x, y, z и в диагональных направлениях координатных плоскостей x = const, y = const, z = const. Методом конечных разностей рассчитаны зависимости от времени перемещений блоков и их скоростей, а также их спектральные плотности при нестационарном нагружении типа «центр расширения». Исследовано влияние вязкости на затухание амплитуды возмущений с ростом времени и расстояния. Проведено сопоставление с решением задачи Шарпа для однородной упругой среды.
Ключевые слова: блочная среда, волновое движение, полость, численное моделирование, вязкоупругое деформирование, дисперсия
Model for the dynamic loading of a cavity in a three-dimensional block medium
N.I. Aleksandrova
Chinakal Institute of Mining SB RAS, Novosibirsk, 630091, Russia
This paper examines wave propagation under dynamic loading of the surface of a cavity in a block medium. A three-dimensional mathematical model of nonstationary viscoelastic deformation of the block medium has been developed. The model is based on the idea that the dynamic behavior of the block medium can be approximately described as the motion of rigid blocks due to the compliance of interlayers between them, and that the interlayer deformation can be approximately described by the Kelvin-Voigt model. Hence the block medium is modeled by a homogeneous three-dimensional lattice of point masses connected by springs and dampers in the x, y, z axis directions and in the diagonal directions of the coordinate planes x = const, y = const, z = const. The finite difference method has been applied to calculate the time dependences of block displacements and velocities, and their spectral densities under nonstationary loading of the "center of expansion" type. The effect of viscosity on wave attenuation with time and distance is studied. The obtained solution is compared to the Sharp problem solution for a homogeneous elastic medium.
Keywords: block medium, wave motion, cavity, numerical modeling, viscoelastic deformation, dispersion
1. Введение
До последнего времени в геомеханике и геофизике широко применялась теория деформирования горного массива как однородной среды, динамика которой описывается хорошо разработанной линейной теорией распространения упругих волн. На основе этой теории построены методы расчета напряженного состояния горных пород вблизи выработок, способы обработки и интерпретации сейсмологической информации в геофизике и горном деле. Серьезный повод к пересмотру сложившихся взглядов дают исследования [1-7], которые свидетельствуют о необходимости учета блочного
строения горных пород в математических моделях. Особенно большое влияние на развитие такого подхода оказала фундаментальная концепция М.А. Садовского [1], согласно которой массив горных пород представляет собой систему вложенных друг в друга блоков разного масштабного уровня, соединенных между собой прослойками, состоящими из более слабых трещиноватых пород.
Как показал анализ размеров блоков в масштабах от фракций породного массива до геоблоков земной коры, отношение размеров соседних по масштабу блоков а = 1Н+Х/ 1Н обладает определенной устойчивостью:
© Александрова Н.И., 2017
а ~ 1.4 [7]. В работе [8], посвященной изучению отношения линейных размеров блоков горных пород к раскрытию трещин между ними в иерархической структуре массивов, экспериментально найден статистический инвариант блочной структуры |Л, равный отношению толщины прослойки между блоками одного масштаба к характерному размеру блока. Например, оценка для сульфидных руд Норильского месторождения [8] дает следующее значение: ц = (0.5-1.0) • 10-2. Исследования деформационных свойств прослоек между блоками разных размеров сейсмическими методами [9] показали, что жесткость прослоек изменяется обратно пропорционально размеру блоков. В [10] показано, что при взрывах разного масштаба линейные размеры консолидированных блоков определяются мощностью взрыва, расстоянием до источника, прочностными характеристиками нарушений сплошности и лежат, как правило, в диапазоне 8-10 м/кт-1/3.
Прослойки между блоками часто представлены более слабыми трещиноватыми породами. Такая структура среды сказывается на процессе распространения волн. Наличие прослоек с ослабленными механическими свойствами приводит к тому, что деформирование блочного горного массива как в статике, так и в динамике происходит в основном за счет их податливости. Как отмечено в работе [6], структура блочной среды является причиной различных динамических явлений, которые отсутствуют в однородной среде, а значит, не могут быть описаны ее моделями. Среди таких динамических явлений выделим распространение маятниковых волн [4, 6], обладающих низкой скоростью распространения (существенно меньшей скорости упругих волн в блоках), большой длиной волны (даже при коротком импульсном воздействии) и слабым затуханием.
Теоретическим аспектам разработки одномерных математических моделей вязкоупругого деформирования блочных сред посвящены работы [ 11-16], в которых показано, что представление блоков как массивных жестких тел позволяет выделить из сложного динамического состояния блочной среды ту ее часть, которая определяется деформированием прослоек между блоками. Такое представление блоков дает возможность достаточно точно описать низкочастотные волны маятникового типа, возникающие при ударном воздействии. В [11-16] показано, что при импульсном воздействии в блочной среде возникает широкополосное возмущение, которое по мере распространения разделяется на высокочастотное, характерное для собственных колебаний блоков, и низкочастотное. Лабораторные эксперименты с использованием физических моделей блочных сред [6, 11-14] показали, что высокочастотные волны быстро затухают и основное сейсмическое воздействие оказывают низкочастотные волны. В отличие от работ [4-14], в которых основное внимание уделено изучению
результатов лабораторных и натурных экспериментов в блочных средах, в работах [15, 16] проведен теоретический анализ влияния блочно-иерархической структуры породного массива на процесс распространения волн и спектральные характеристики возмущений.
Двумерные модели блочных сред рассматривались в работах [17-35]. В [17, 18] блочная среда моделируется как упаковка из круглых относительно прочных дисков (клеточных автоматов) одинакового радиуса, разделенных границами со свойствами, существенно отличающимися от свойств блоков. В [19, 20] блочная среда моделируется набором шестиугольных элементов, взаимодействующих с силой равной отношению давления к площади соприкосновения. Для данной модели была изучена динамика ненапряженного полупространства при условии, что слой поверхностных элементов подвержен гармоническим колебаниям.
В [21] исследовались некоторые особенности сейсмического процесса в блочной среде на модели, описывающей нерегулярную упаковку жестких дисков, радиусы которых распределены в интервале 0.65-1.35. В качестве меры сейсмической активности использовалась кинетическая энергия блоков. Исследовалась зависимость суммарной кинетической энергии блоков от времени нагружения. Изучалось также влияние микроколебаний на процесс выделения энергии при деформировании блочной среды. В [22, 23] горный массив представлен блоками кубической формы, расположенными по одну сторону от трещины. Причем больший структурный блок может состоять из меньших блоков. Между собой и массивом блоки взаимодействуют по закону сухого трения. Для такой модели оценены параметры движения блоков при действии плоской волны давления треугольной формы.
В работах [24-26] исследовалось динамическое поведение двумерной блочной среды. В этих работах был использован такой же подход, который применялся в [15] для одномерных блочных сред. В [24-26] предполагалось, что блоки прямоугольной формы взаимодействуют через податливые прослойки. В [24] блоки предполагались жесткими, в [25, 26] — упругими. В [26] для моделирования двумерной блочной среды использовались также уравнения ортотропного континуума Кос-сера.
Более простую модель двумерной блочной среды можно получить, если считать блоки сосредоточенными массами, соединенными пружинами и демпферами. В этом случае регулярную блочную среду можно представить как периодическую решетку масс. Различные варианты такой модели использованы для описания плоского и антиплоского деформирования дискретных сред в работах [27-35]. Впрочем, в работе [35] основное внимание уделяется непрерывным моделям, описывающим дискретные микроструктуры.
Несмотря на очевидную ограниченность возможностей использования периодических решеток как моделей для описания реальных блочных сред, их преимуществом является возможность привлекать аналитические и численные методы, а также качественно описывать динамические явления, свойственные этим средам. Изначально теория волн в периодических структурах возникла в динамической теории кристаллических решеток, затем она нашла применение в механике композитов, но долгое время не использовалась для описания динамических событий в горном массиве. Однако исследование распространения волн в периодических блочных структурах интересно также для приложений в сейсмике.
В данной работе блочная среда моделируется как трехмерная решетка масс, соединенных пружинами и демпферами в осевых и диагональных направлениях координатных плоскостей. Ниже получено конечно-разностное решение задачи для решетки при нестационарном воздействии типа «центр расширения» и проведено сопоставление с аналитическим решением задачи для сферической полости в однородной упругой среде [36].
2. Постановка задачи и уравнения движения
Исследуется нестационарная пространственная задача о воздействии типа «центр расширения» в блочном массиве. Блочная среда моделируется однородной трехмерной решеткой, состоящей из точечных масс, соединенных пружинами и демпферами в направлениях осей x, y, z и в диагональных направлениях плоскостей x = = const, y = const, z = const, как показано на рис. 1, a. Здесь использованы обозначения: u, v, w — перемещения в направлениях осей x, y, z; n, m, k—номера блоков в направлениях x, y, z. Пусть длина пружин в направлениях осей x, y, z равна l. Положим l = 1. Пусть точка O соответствует началу координат (рис. 1, б) и в точках O (0,0,0), Ai (1,0,0), A2 (1,1,0), A3 (0,1,0), A4 (1,0,1), A5 (1,1,1), A6 (0,1,1), A7 (0, 0,1) приложены силы, равные по величине и направленные от центра симметрии
(точка О1) куба с вершинами в точках О, А1, А2, А3, А4, А5, А6, А7, как показано на рис. 1, б. Пусть, кроме того, в решетке отсутствуют связи между массами, находящимися в вершинах этого куба. На рис. 1, б эти связи отмечены штриховыми линиями. Такое нагружение моделирует воздействие типа «центр расширения», приложенное к полости в блочной среде. Приближенно можно считать, что эта полость имеет радиус а = л/3/ 2 (радиус сферы, описанной вокруг единичного куба).
Для теоретического описания деформирования пружин и демпферов используется реологическая модель Кельвина-Фойгта [37]. С учетом обозначений
Л пп/п,т,к = Л+1, - 2/п, т,к + /п-1,т,к,
Фпт/п,т,к = /п+1,т+1,к + /п-1,т-1,к - 4/п,т,к +
+ /п+1,т-1,к + /п-1,т+1,к, (1)
^ пт/п, т, к = /п+1, т+1, к + /п-1, т-1, к /п+1, т-1,к
fn-1,;
m+1, k,
уравнения движения массы с координатами п, т, к имеют вид [38]:
Мйп
,т ,к к1Лппип ,т ,к + к2[(Ф пк + Ф пт ^ип ,т,к + + ^птУп,т,к + ^пк™п,т,к V2 + Х1Лппип,т,к + + Х2[(Фпк + Фпт )ип,т,к пт^п,т,к +
+ ^пк^п,т,к V2,
Мп,т ,к = к1Л ттУп ,т,к + к2[^ птип ,т,к +
+ (Ф 7+ф )у 7 +Т ,1/2 +
\ тк пт' п,т,к тк п,т,кJ/
• • (2)
+ (Фтк + Ф
пт ,т,к + * тк^п,т,к V2, М^п,т,к = к1Л кк^птк + к2[^ пкип,т,к + + (Фтк + Фпк ) ™п,т,к V2 + + Х1Лкк^п,т,к + Х2[^пкип,т,к + ^тЛ,т,к + + (Ф тк + Ф пк ) ™п ,тк У2'
Здесь к1, к2, Л15 X2 — жесткости пружин и вязкости демпферов в осевых и диагональных направлениях соответственно. Параметры М, к1, к2, Х1, X2, входя-
Рис. 1. Схема соединения масс пружинами и демпферами в трехмерной модели блочной среды (а); схема приложения нагрузки (б)
щие в уравнения (2), имеют одни и те же значения во всех точках среды. Начальные условия для уравнений (2) нулевые:
t = 0: ип
= и1п,т,к = 0,
= Уп т к = Щ т к = Щ
= 0.
(3)
Уп,т,к Уп,тк Щп,тк Щп,тк В частном случае, когда Х1 = X 2 = 0 , модель блочной среды (2) совпадает с частным случаем модели решетки Борна-Кармана [39], в которой не учитывается сдвиговая жесткость.
Поскольку плоскости х = 0.5, у = 0.5 и г = 0.5 являются плоскостями симметрии волнового процесса, то на них выполняются следующие условия:
Чт,к = и1,т,к, У0,т,к = У1,т,к, Щ,т,к = Щт,к, ип,0,к = ип,1,к, Уп,0,к = -уп,1,к, Щп,0,к = Щп,1,к,
(4)
ип,т,0 ип,т1, Уп,т,0 Уп,т,1, Щп,т,0
*п,т,1.
Далее, в силу симметрии, будем рассматривать процесс распространения волн только в области: п > 0, т >0, к > 0. Пусть в точке ^(1,1,1) (рис. 1, б) приложена нагрузка О = (6х ,6у, 6г) = (1,1,1)6 (), где 6 (t — ее амплитуда в момент времени t.
В (2) заменим каждое разностное соотношение (1) его дифференциальным приближением, например
Л ппип,т,к = ип+1,т,к 2ип,т,к + и г4 э4
п-1,т,к
Д2 >4 Д4
~ 12 " и + I Э и +
~ эХ2+йдХ4+....
Полагая I ^ 0, опустим слагаемые порядка выше чем 12. В результате получим, что при переходе к сплошной среде уравнениям (2) при Х1 = Х2 = 0 соответствуют уравнения изотропной однородной упругой среды:
Ми = I
„, .Э2и (к1 + 2к2^-т + к2 Эх
С Э2и Э2и ) Эу2 + Эг 2
С
+2к2
Э2у Э2щ +
МУ = I
ЭхЭу ЭхЭг
\
.. ЧЭ 2у
(к1 + 2к2)—Т + к2
Эу2
С
+2к
Э 2и Э 2щ +
МЩ = I
ЭхЭу ЭуЭг
д^щ ЭГ
(к1 + 2к2)ГТХ + к2
Э 2у Э" у Эх2 Эг2
С Э 2 Щ д^Щ ) 2 + Эг2
(5)
Эу
+2к
Сд 2
и Э 2 у ^
ЭхЭг ЭуЭг
\
с коэффициентом Пуассона а = 0.25. Это соответствует следующим скоростям продольных Ср и сдвиговых с8
Ср = I,
к1 + 2к2 С = , ¡К
М 'С Чм •
3. Дисперсионные свойства
Для того чтобы получить дисперсионные соотношения модели, описываемой уравнениями (2), применим к каждому из этих уравнений преобразование Лапласа по времени t (параметр р, значок Ц) и дискретные преобразования Фурье по переменным п, т, к (параметры qx, qy, qz соответственно и значки Fn, Fm, Fk). Введем обозначения:
\LFFF, Тг ,,
и = (ия. Ж = Щ
)
V = (Уп т к ) М=
Система преобразованных по Лапласу и Фурье уравнений (2) в случае Х1 = Х2 = 0 имеет вид:
С
а11 + М
а21
а31
а12 «22 + ММ
а32
а13
а23
а33
+ М
и \
V Ж
С 0 ^
0
V У
где
волн:
а11 = 42 + 4Y + 4X(К + 2 - 22 - 2Y), а22 = 42 + 4X + 4Y(К + 2 - 22 - 2X), а33 = 4Y + 4 X + 42 (К + 2 - 2Y - 2 X), а122 = а2 = 64 X (1 - X )Y (1 - Y), а2 = а321 = 64X (1 - X )2(1 - 2), а|3 = а322 = 64Y (1 - Y)2 (1 - 2), X = ^т2(Ях1!2), Y = ¿т2^уЦ2), 2 = sin2(qz^2), к^к2 = К.
Делая замену р = в дисперсионном операторе А = det(А + ММи), где А — матрица, элементы ау которой определены выше, J — единичная матрица, и приравнивая его нулю, получим дисперсионное уравнение для определения зависимости частоты ю от волновых чисел дх, Цу, :
А = det( А-М ю2/Д2) =0. (7)
Поскольку матрица А симметрична, то все ее соб-ственн^1е числа вещественны. Численные расчеты показывают, что все главные миноры матрицы А неотрицательны при Х1 = Х2 = 0. Поэтому все корни дисперсионного уравнения (7) также неотрицательны, т.е. Мю2/к2 > 0. Следовательно, частота ю имеет только вещественные значения, если Х1 =Х 2 = 0.
Векторы фазовых и групповых скоростей и их модули определяются формулами:
^ ^^-^-2 (Цх. Чу , Чг )>
а ^ + д „ + а.
ш
I с, I =
I 2 . 2 . 2 '
^Чх + Чу + ^
Эш Эш Эш
(8)
Эш ^ 2 Эш ^ 2 г Эш Л
с А = А + +
ё1 А 1Эч*) дЧу \ у ) ^ J
Пусть = 0. Тогда 7 = 0, а
13
23
- 0 и определи-
тель (7) распадается на два сомножителя: Д = (а33 -Мш2/к2)[(а11 -Мш2/k2)х
х (а22 - Мш2/k2) - а22 ] = 0. Это уравнение имеет следующие корни: ш12 = ш0 [2( X + Y)(K + 3) - 8 ^ ±
± 2{[(X + Y)(K +1) - 4X7]2 -
- 4 Х7 [(К +1 - 27)(К +1 - 2 X ) -
- 4(1 - 7 )(1 - X )]}1/2]1/2,
(9)
ш
= 2ш0л/ X + 7,
где ш0 .
Обозначим через сг j вектор фазовой скорости, отвечающий частоте Шj, где j = 1, 2, 3. Используя формулы (8), (9), можно вычислить значения | сг j | при Ч2 = 0. В первую очередь нас интересуют скорости длинных волн. Поэтому приведем значения | сг j |, полученные в предположении, что чх, чу ^ 0 и д2 = 0:
Чх = 0 Чу ^ 0:
| СГ,1
= ср =1
k1 + 2к2
М
^,2
= | с
= Сс 1 = с
^3 1 _ с8,1
(10)
Чх = Чу ^ 0:
|си|
| Сг,2 |
\сы
Модули фазовых и групповых скоростей для произвольных значений Чх, Чу, Ч2 определялись численно из уравнения (7) с помощью формул (8). Эти расчеты показали, что модули фазовых и групповых скоростей длинных волн в блочной среде совпадают. Из этого факта следует, что бесконечно длинные волны распространяются в блочной среде без дисперсии и формируют низкочастотные маятниковые волны.
На рис. 2 представлены графики зависимости модуля фазовой скорости от волновых чисел Чх, Чу в сечении ч2 = 0 при различных соотношениях жесткости пружин в осевых и диагональных направлениях. Расче-
Рис. 2. Зависимость модуля фазовой скорости от волновых
чисел
при ч2 = 0
ты выполнены для решетки с параметрами М = I = 1, X = X2 = 0. На каждом графике (рис. 2) верхняя поверхность соответствует фазовой скорости продольных волн, нижние две поверхности — фазовым скоростям сдвиговых волн. В формулах (10) и на рис. 2 видно, что при различных соотношениях волновых чисел модули фазовых скоростей бесконечно длинных волн в трехмерной блочной среде могут быть разными, если к1 ф к2.
Как видно из формул (10), скорость бесконечно длинных продольных и сдвиговых волн в блочной среде зависит от длины пружин I, массы блоков Ми жесткос-тей к1, к2 осевых и диагональных пружин. Физически, применительно к массиву, параметр «длина пружины I» соответствует размеру блоков. Это показано в [15] на примере одномерной цепочки упругих стержней, соединенных пружинами. Если в уравнениях (2) положить жесткость диагональных связей равной нулю (к2 = 0), то полученные уравнения будут описывать распространение одномерных продольных волн в направлениях осей х, у, г. При этом сдвиговые волны формироваться не будут.
В [38] показано, что для решетки, представленной на рис. 1, а, существуют только три собственных значения частоты коротковолновых возмущений в случае к1 = к2 и Х1 = X 2 = 0:
2ю0,2л/2ю0,2^/2ю0. (11)
Расчеты показали, что модуль групповой скорости, соответствующий этим значениям, равен нулю, и, следовательно, можно ожидать, что эти значения частот будут резонансными для данной модели блочной среды.
Как следует из дисперсионных соотношений (9) и расчетов, представленных на рис. 2, первыми движутся низкочастотные продольные волны, имеющие максимальную скорость в блочной среде. Позади них движутся низкочастотные сдвиговые волны и коротковолновые возмущения.
4. Устойчивость разностной схемы
Уравнения (2) с условиями (3) и (4) решались методом конечных разностей по явной схеме. Для вторых производных по времени использовалась центрально-разностная аппроксимация второго порядка точности:
^+1
^-1
ип,т,к :
ип,т,к 2ип,т,к + ип,т,к
s = 0, 1, 2,
Для первых производных по времени использовалась разность «назад» первого порядка точности:
5 _ 5-1
ип,т,к ип,т,к п 1
= —т-—, s = 0, 1, 2, ....
и
т
Здесь т — шаг разностной сетки по времени; ихптк — значение перемещения ип тк (}) в момент времени t = = лт, s — номер слоя по времени в конечно-разностной схеме. Для перемещений уп-,к ^) , щпт к (t) использовались аналогичные аппроксимации. Можно показать, что условие устойчивости разностных уравнений, соответствующих уравнениям (2), имеет вид:
С I—г;;-^
т < min
Х1 + Х2 к1 + к2
Х
-1 +
1+Мк2
2Х2
-1+. 1 +
М (к1 + к2)
(Х1 +Х2)2
Х1 к1
-1 + 1 +
Мк1
Х12
Х1 + 2Х:
С
к1 + 2к2
-1 + . 1 +
М (к1 + 2к2)
м
(Х1+2Х2)2
\ / ^
Конечно-разностные расчеты подтвердили справедливость этого условия.
5. Результаты численных расчетов. Ступенчатая нагрузка
В этом разделе представлены результаты численных расчетов возмущений при действии ступенчатой нагрузки 6 (ь) = 60 Н^), приложенной, как показано на рис. 1, б. Здесь И(Ь) — функция Хевисайда; 60 — амплитуда нагрузки.
Для сравнительного анализа результатов, полученных по модели блочной среды (2) в случае ступенчатого воздействия типа «центр расширения», нам потребуется решение динамической задачи о расширении сферической полости в упругой среде под действием ступенчатой радиальной нагрузки, приложенной к поверхности полости.
Пусть в неограниченном упругом однородном и изотропном пространстве имеется полость радиуса а0. Пусть на ее границе заданы радиальные напряжения постоянной амплитуды: аг = -р0 (ь > 0, г = а0). В сферических координатах уравнение движения упругой среды относительно радиальных перемещений и имеет вид:
Эt
2 = с2
Ср2 =-
С Э
эГ 2
Е(1 -а)
Э и 2Эи „и
г Эг г
р(1 + а)(1 - 2 а)
Здесь г — радиальная координата; Ср — скорость продольных волн в упругой среде; Е — модуль Юнга; р — плотность упругой среды; а — коэффициента Пуассона. Начальные условия нулевые.
Решение данной задачи впервые получено в связи с расчетом параметров взрывных волн Дж. Шарпом [36] в случае а = 0.25. Для произвольных значений коэффициента Пуассона оно имеет вид:
и(г, t) =
Р0 а0(1 + а) I а0
2гЕ
- 2 IV1- 2а sin (со£) +
г
Н 0(^),
(12)
и (г, г) =
г рю га 1
С1 га^
а0 Vl - 2а
cos(сb £)
sin (оо £) +
Н 0©,
т
Рис. 3. Зависимости перемещений от времени. Кривые 1 — решение (12) для полости в упругой среде, для полости в блочной среде: кривые 2 — Х1 = Х 2 = 0, кривые 3 — Х1 = X 2 = 0.05, кривые 4 — Х1 =Х 2 = 0.1
где
а = -
1 - 2а 1 -а '
со = -
У1 - 2а Ср 1 -а а0
Р = аСЕ-.
Определим жесткость пружин, при которой скорости продольных волн в блочной (случай к1 = к2) и упругой средах совпадают: Е (1 -а)
, (13)
(р(1 + а)(1 - 2а)
При I ^ 0 уравнения (2) переходят в уравнения (5), описывающие движение однородной изотропной упругой среды с коэффициентом Пуассона а = 0.25. Поэтому положим в (13) а = 0.25. Пусть р = 1 и Е = 1, М = 1, I =1. В результате из (13) следует, что к1 = 0.4.
В [40] «кажущийся» радиус а0 и «кажущееся» давление р0 упругой среды в теоретическом решении (12) подбирались из условия совпадения максимумов радиальных перемещений в формулах (12) и в экспериментальных данных для каждой точки по радиусу, в которой проводилось сравнение. В отличие от [40], в данной работе мы подбирали «кажущееся» давление р0 в упругой модели из условия совпадения остаточных радиальных перемещений при больших значениях вре-
мени для блочной и упругой сред, а радиус полости а0 полагали равным л/э/2. В результате было получено, что р0 - 2.177 Q0. Это значение р0 является общим для каждой точки по радиусу в упругой модели. Пусть г — расстояние от центра симметрии блочной среды, имеющего координаты (0.5, 0.5, 0.5), до точки с координатами (п, т, к) (рис. 1, б):
г = у1(п - 0.5)2 + (т - 0.5)2 + (к - 0.5)2 . Это значение радиуса использовалось в формулах (12) для определения перемещений, скоростей и ускорений упругой среды.
Все расчеты, представленные ниже, получены при следующих значениях параметров блочной среды: Q0 = = М = I = 1, к1 = к2 = 0.4, и упругой среды: р0 = 2.177, р = 1, Е = 1, а = 0.25. Поскольку волновой процесс симметричен относительно плоскостей х = 0.5, у = 0.5,
2 = 0.5, то конечно-разностные расчеты для блочной среды проводились только для области п > 1, т > 1, с учетом условий (4) и с шагом по времени т = 0.25.
На рис. 3, 4 приведены результаты расчетов зависимости перемещений ип11 и скоростей перемещений ип1д от времени для блочной среды при различных значениях параметра вязкости (кривые 2 — Х1 = X 2 = 0, кривые 3 — Х1 = X 2 = 0.05, кривые 4 — Х1 =Х 2 = 0.1) и различных значениях координаты п (п = 16, 46). На
Рис. 4. Зависимости от времени скорости перемещений ип11. Х1 =Х2 = 0 (тонкие линии), 0.05 (линии средней толщины), 0.1 (толстые линии)
рис. 3 кривые 1 соответствуют решению (12) для полости в упругой среде. Вертикальные линии на рис. 3, 4 соответствуют моментам времени прихода фронтов продольных и сдвиговых волн в точку с координатами Оп, 1, 1): гр — (п-1)/Ср, г8 = (п- 1)/с8, где ср, с определены в (6).
На рис. 3 видно, что в решении для блочной среды появляются высокочастотные осцилляции, отсутствующие в решении (12) для упругой среды, и решение для блочной среды существенно отличается от решения для упругой среды. Тем не менее есть качественное соответствие, а именно: перемещение блочной среды осциллирует относительно перемещения упругой среды.
Рисунки 3, 4 показывают, что при "1 — К2 — 0 за фронтом сдвиговой волны появляются высокочастотные осцилляции, отношения максимальных амплитуд которых к максимальным амплитудам низкочастотных волн увеличиваются с удалением от полости, при этом абсолютные значения максимумов их амплитуд падают. Этот же вывод подтверждают графики спектральной плотности скоростей перемещений т
G (&, ю) — / и (0е_йог dt о
в различных точках по пространству (п = 16, 31, 46) для различных значений параметра вязкости, представленные на рис. 5. Остальные параметры те же, что на рис. 3, 4. Спектральная плотность вычислена с помощью быстрого преобразования Фурье (Т = 150). На рис. 5 вертикальные линии соответствуют собственным частотам (11) коротковолновых возмущений в блочной среде. Видно, что на близких расстояниях от места воздействия есть пики на графиках спектров скоростей в точках, соответствующих частотам (11). С удалением от места воздействия эти пики исчезают и остаются две несущих частоты: ю « 1 и ю « 2. В случае упругого деформирования блочной среды (т.е. К1 — К 2 — 0) с ростом координаты п происходит перестройка волнового процесса: вклад колебаний с частотой ю « 1 уменьшается и основными становятся колебания с частотой ю « 2. В случае вязкоупругого деформирования высокочастот-
ные колебания с несущей частотой ю « 2 очень быстро затухают с ростом координаты п и при больших значениях п остаются только колебания с частотами ю< 2^М -1.265 и несущей частотой ю-
2 М/ (5^) — 1. С ростом параметра вязкости максимальная амплитуда спектральной плотности в окрестности ю « 1 также быстро затухает.
Таким образом, как показывают рис. 3-5, наличие диссипативных свойств у прослоек приводит к дополнительному затуханию волн с частотой ю « 1 и к очень быстрому затуханию высокочастотных волн с частотой ю « 2. На рис. 3, 4 видно также, что в случае К1 Ф 0, К2 Ф 0 низкочастотные маятниковые волны, распространяющиеся с максимальной скоростью в блочной среде, затухают существенно медленнее, чем высокочастотные колебания. Это приводит к тому, что в блочных горных массивах на больших расстояниях от места воздействия основной вклад в волновой процесс вносят низкочастотные маятниковые волны.
Были проведены расчеты максимальных амплитуд осциллограмм возмущений и, и, и от координаты п для различных значений параметра вязкости: К1 — 0.00, 0.01, 0.05, 0.07, 0.10. Для каждой кривой была найдена ее аппроксимация степенной функцией вида / — Апа, где константы А и а подбирались так, чтобы значения функции / в точках п = 5, 10, 15, ..., 50, 55, 60 были наиболее близки к значениям соответствующих функций и, и, и в тех же точках. Для констант а была найдена их аппроксимация линейной функцией вида а — -(аК1 + Ь). Анализ численных результатов показал, что максимальные амплитуды и, и, и как функции переменной п, стремятся к нулю с ростом п. Более того, скорость убывания возрастает с ростом вязкости прослоек и определяется формулами:
тах|ипП(г)| «А(ц)п (3'46К1+ г
тах |ипП(г)| « А(&)п~(7'т'+1'2),
г
тах |1 1(г)|« А(й)п-(11'32"+1'2) .
г ''
Рис. 5. Зависимости спектральной плотности скорости перемещений от частоты. К1 — К2 — 0 (тонкие линии), 0.01 (линии средней толщины), 0.05 (толстые линии)
На рис. 6 представлены графики распределения радиальной скорости перемещений
. Vпт,1(т - а5) + ипт,1(п - 0.5)
иг = , -
д/(п - 0.5)2 + (т - 0.5)2 + 0.52 как функции переменных п, т в сечении k = 1 в момент времени t = 100 при действии ступенчатой нагрузки. Расчеты выполнены для решетки с параметрами Q0 = = М = 1=1, Х1 =Х2 = 0. Жесткость диагональных и осевых пружин варьировалась. Видно, что если к1 = к2,
к1 = 0.8, к2 = 0.2
то фронт продольной волны является окружностью, т.е. при к1 = к2 имеем «изотропную» решетку: низкочастотные продольные волны распространяются с одинаковой скоростью во всех направлениях. Если к1 > к2, то скорость распространения возмущений вдоль осей больше, чем вдоль диагоналей, а если к1 < к2, то скорость распространения возмущений вдоль диагоналей больше, чем вдоль осей. Видно, что если жесткости пружин в осевых и диагональных направлениях отличаются значительно, то и фронт продольной волны значительно отличается от окружности. Анализ амплитуд возмущений показал, что если к1 > к2, то амплитуда радиальных скоростей на оси больше, чем на диагонали; если к1 < к2, то наоборот.
6. Результаты численных расчетов. Импульсная и монохроматическая нагрузки
В этом разделе представлены результаты численных расчетов перемещений в задаче Шарпа для блочной среды при импульсном нагружении длительности t0, моделируемом полусинусоидой:
Q^) = QoSin(шt)Н ^)Н (^ -1), ш = П(14) и при воздействии монохроматической нагрузки с частотой ш:
Q (t) = QoSin(шt) Н ^). (15)
Нагрузки (14), (15) приложены, как показано на рис. 1, б.
На рис. 7 представлены результаты расчетов зависимости перемещений ип11 от времени для блочной среды. Тонкие линии соответствуют нагрузке (15), толстые линии — нагрузке (14). Вертикальные пунктирные линии соответствуют времени прихода фронтов продольных и сдвиговых волн в точку с координатами (п, 1, 1). Все расчеты проведены для двух значений ш = 2 и 1 и следующих значений остальных параметров: £0 = М = I = 1, к1 = к2 = 0.4, Х1 =Х 2 = 0, п = 21, т = 0.25.
Рис. 6. Графики распределения радиальной скорости перемещений по переменным х, у в сечении k = 1
Рис. 7. Зависимости перемещений от времени при действии импульсной (толстые линии) и монохроматической (тонкие линии) нагрузок в случае ш = 1 (а) и 2 (б)
Видно, что амплитуда колебаний при монохроматическом воздействии существенно выше, чем при импульсном воздействии. Тем не менее резонансного роста возмущений, который можно было бы ожидать на основе анализа спектральных характеристик, не наблюдается. Проводились также расчеты при монохроматическом воздействии с собственной частотой коротковолновых возмущений, и также явление резонанса не наблюдалось. Заметим, что в случаях одномерных и двумерных блочных сред при воздействии с частотой собственных колебаний резонансный рост возмущений всегда был. По-видимому, в трехмерном случае пространственное рассеяние настолько велико, что даже при монохроматическом воздействии с собственной частотой резонанс не может реализоваться.
7. Заключение
Предложена трехмерная математическая модель блочной горной породы. Эта модель основана на представлении, что динамическое поведение блочной среды может быть приближенно описано как движение жестких блоков за счет податливости прослоек между ними и что деформацию прослоек можно приближенно описать моделью Кельвина-Фойгта.
С использованием данной модели численно решена задача о нестационарном воздействии типа «центр расширения» на поверхность полости в блочной среде. В результате показано, что основной вклад в волновой процесс в блочной среде вносят низкочастотные маятниковые волны в районе фронта продольной волны и за ними движутся коротковолновые возмущения с частотой близкой к частоте поверхностных волн блочной среды.
Показано, что наличие блочной структуры в среде приводит к следующим изменениям поведения блочной среды по сравнению с тем, что предсказывает модель однородной упругой среды, получаемая усреднением механических свойств блочной среды:
- в блочной среде распространяются низкочастотные слабо затухающие маятниковые волны со скоростями, которые определяются массой блоков, их размерами и зависят от реологических свойств прослоек, за ними движется группа высокочастотных колебаний;
- волны в блочной среде распространяются с дисперсией, которая отсутствовала в однородной среде;
- наличие диссипативных свойств у прослоек приводит к быстрому затуханию высокочастотных волн и к дополнительному затуханию низкочастотных волн.
Эти результаты свидетельствуют о необходимости учета блочной структуры горных пород и реологических свойств прослоек при расчете сейсмических волн.
Работа выполнена при финансовой поддержке проекта ОНЗ РАН-3.1.
Литература
1. Садовский М.А. Естественная кусковатость горной породы // ДАН СССР. - 1979. - Т. 247. - № 4. - С. 829-832.
2. Садовский М.А., Болховитинов Л.Г., Писаренко Б.Ф. О свойстве дискретности горных пород // Изв. АН СССР. Физика Земли. -1982. - № 12. - С. 3-18.
3. Садовский М.А., Писаренко Б.Ф. Сейсмический процесс в блоковой среде. - М.: Наука, 1991. - 96 с.
4. Курленя М.Б., Опарин Б.Н., Бостриков Б.И. О формировании упругих волновых пакетов при импульсном возбуждении блочных сред. Волны маятникового типа U^// ДАН СССР. - 1993. - Т. 333. -№ 4. - С. 3-13.
5. Курленя М.Б., Опарин Б.Н., Бостриков Б.И. Волны маятникового типа. Ч. I: Состояние вопроса и измерительно-вычислительный комплекс // ФТПРПИ. - 1996. - № 3. - С. 3-7.
6. Курленя М.Б., Опарин Б.Н., Бостриков Б.И. Волны маятникового типа. Ч. II: Методика экспериментов и основные результаты физического моделирования // ФТПРПИ. - 1996. - № 4. - С. 3-38.
7. Курленя М.Б., Опарин Б.Н., Бостриков Б.И., Аршавский Б.Б., Ма-мадалиев Н. Волны маятникового типа. Ч. III: Данные натурных измерений // ФТПРПИ. - 1996. - № 5. - С. 3-27.
8. Курленя М.Б., Опарин Б.Н., Еременко А.А. Об отношении линейных размеров блоков горных пород к величинам раскрытия трещин в структурной иерархии массива // ФТПРПИ. - 1993. - № 3. -С. 3-21.
9. Костюченко Б.Н., Кочарян Г.Г., Павлов Д.Б. Деформационные характеристики межблоковых промежутков различного масштаба // Физ. мезомех. - 2002. - Т. 5. - № 5. - С. 23-42.
10. Кочарян Г.Г., Спивак А.А. Движение блоков горной породы при крупномасштабных подземных взрывах. Ч. 1: Экспериментальные данные // ФТПРПИ. - 2001. - № 1. - С. 71-83.
11. Александрова Н.И., Шер Е.Н. Моделирование процесса распространения волн в блочных средах // ФТПРПИ. - 2004. - № 6. -С. 49-57.
12. Александрова Н.И., Черников А.Г., Шер Е.Н. Экспериментальная проверка одномерной расчетной модели распространения волн в блочной среде // ФТПРПИ. - 2005. - № 3. - С. 46-55.
13. Александрова Н.И., Черников А.Г., Шер Е.Н. О затухании маятниковых волн в блочном массиве горных пород // ФТПРПИ. - 2006. -№ 5. - С. 67-74.
14. Александрова Н.И., Шер Е.Н., Черников А.Г. Влияние вязкости прослоек на распространение низкочастотных маятниковых волн в блочных иерархических средах // ФТПРПИ. - 2008. - № 3. -С. 3-13.
15. Александрова Н.И. О распространении упругих волн в блочной среде при импульсном нагружении // ФТПРПИ. - 2003. - № 6. -С. 38-47.
16. Шер Е.Н., Александрова Н.И., Айзенберг-Степаненко М.Б., Черников А.Г. Влияние иерархической структуры блочных горных пород на особенности распространения сейсмических волн // ФТПРПИ. - 2007. - № 6. - С. 20-27.
17. Астафуров С.Б., Шилько Е.Б., Псахье С.Г. Влияние стесненных условий на характер деформирования и разрушения блочных сред при сдвиговом нагружении // Физ. мезомех. - 2009. - Т. 12. - № 6.-С. 23-32.
18. Псахье С.Г., Смолин А.Ю., Коростелев С.Ю., Дмитриев А.И., Шилько Е.Б., Астафуров С.Б. Метод подвижных клеточных автоматов и его применение при моделировании на разных масштабах // Механика — от дискретного к сплошному / Под ред. В.Ф. Фомина. - Новосибирск: Изд-во СО РАН, 2008. - С. 88-128.
19. Бенгрович Д.Б., Губар И.М., Шеремет Г.П. Распространение нелинейных волн в дискретной упругопластической среде //
Материалы XXIV Межд. науч. школы им. академика С.А. Хрис-тиановича «Деформирование и разрушение материалов с дефектами и динамические явления в горных породах и выработках», 22-28 сентября 2014 г., Алушта. - 2014. - С. 46-50.
20. Венгрович Д.Б. Сейсмические волны в блочных средах // Динамика сплошной среды: Сб. научн. тр. - Новосибирск: ИГиЛ СО РАН, 1995. - Вып. 110. Акустика неоднородных сред.- С. 57-62.
21. Кочарян Г.Г., Федоров А.Е. Об особенностях механики сейсмического процесса в блочной геофизической среде // ДАН СССР. -1990. - Т. 315. - № 6. - С. 1345-1349.
22. Кочарян Г.Г. Модель необратимого деформирования горного массива блочной структуры при взрывном воздействии // Взрывное дело. - М.: Недра, 1990. - № 90/47. - С. 30-42.
23. Кочарян Г.Г., Спивак А.А., Будков А.М. Движение блоков горной породы при крупномасштабных подземных взрывах. Ч. 2. Оценки по аналитическим моделям, численные расчеты и сравнительный анализ теоретических и экспериментальных данных // ФТПРПИ. -2001. - № 2. - C. 37-57.
24. Сарайкин В.А. Расчет волн, распространяющихся в сборке из прямоугольных блоков // ФТПРПИ. - 2008. - № 4. - С. 32-42.
25. Сарайкин В.А. Учет упругих свойств блоков в низкочастотной составляющей волны возмущений, распространяющейся в двумерной среде // ФТПРПИ. - 2009. - № 3. - С. 9-24.
26. Садовский В.М., Садовская О.В., Варыгина М.П. Анализ резонансного возбуждения блочной среды на основе уравнений мо-ментного континуума Коссера // РЭНСИТ. - 2013. - Т. 5. - № 1.-С. 111-118. - http://en.rensit.ru/vypuski/116/.
27. Ayzenberg-Stepanenko M. V., Slepyan L.I. Resonant-frequency primitive waveforms and star waves in lattices // J. Sound Vib. - 2008. -V. 313. - P. 812-821. - doi 10.1016/jjsv.2007.11.047.
28. Александрова Н.И., Айзенберг-Степаненко М.В., Шер Е.Н. Моделирование распространения упругих волн в блочной среде при импульсном нагружении // ФТПРПИ. - 2009. - № 5. - С. 2233.
29. Александрова Н.И. Асимптотическое решение антиплоской задачи для двумерной решетки // ДАН. -2014. - Т. 455. - № 1. - C. 3437.
30. Александрова Н.И., Шер Е.Н. Распространение волн в двумерной периодической модели блочной среды. 4.1: Особенности волнового поля при действии импульсного источника // ФТПРПИ. -2010. - № 6. - С. 60-72.
31. Aleksandrova N.I. The discrete Lamb problem: Elastic lattice waves in a block medium // Wave Motion - 2014. - V 51. - P. 818-832. -doi 10.1016/j.wavemoti.2014.02.002.
32. Александрова Н.И. Плоская задача Лэмба для блочной среды с учетом вязкости // Вестник Инженерной школы ДВФУ. - 2014. -№ 3. - C. 39-45. - http://vestnikis.dvfu.ru/vestnik/2014/3/4/.
33. Alessandrini B., Raganelli V. The propagation of elastic waves in a discrete medium // Eur. J. Mech. Solids. A. -1989. - V. 8. - No. 2. -P. 129-160.
34. Jensen J.S. Phononic band gaps and vibrations in one- and two-dimensional mass-spring structures // J. Sound Vib. - 2003. - V 266. -P. 1053-1079. - doi 10.1016/S0022-460X(02)01629-2.
35. Andrianov I.V., Awrejcewicz J., Weichert D. Improved continuous models for discrete media // Math. Probl. Eng. - 2010. - 35 p. - doi 10.1155/2010/986242.
36. Sharpe J.A. The production of elastic waves by explosion pressure // Geophysics. - 1942. - V.7. - No. 2. - P. 144-154.
37. Bland D.R. The Theory of Linear Viscoelasticity. - Oxford: Pergamon Press, 1960.
38. Aleksandrova N.I. Seismic waves in a three-dimensional block medium // Proc. R. Soc. A. - 2016. - V. 472. - No. 2192. - doi 10.1098/ rspa.2016.0111.
39. BornM., von Karman Th. Uber Schwingungen in Raumgittern // Phys. Z. - 1912. - V. 13. - P. 297-309.
40. Healy J.H., Chi-Yu King, O'Neill M.E. Source parameters of the salmon and sterling nuclear explosions from seismic measurements // J. Geophys. Res. - 1971. - V. 76. - No. 14. - P. 3344-3355.
Поступила в редакцию 01.07.2016 г., после переработки 31.08.2016 г.
Сведения об авторе
Александрова Надежда Ивановна, д.ф.-м.н., доц., внс ИГД СО РАН, nialex@misd.ru
