Численное моделирование процесса накопления межблоковых перемещений при низкоамплитудных динамических воздействиях Текст научной статьи по специальности «Физика»

УДК 550.34.01, 550.348

Численное моделирование процесса накопления межблоковых перемещений при низкоамплитудных динамических воздействиях

А.М. Будков, Г.Г. Кочарян, Д.В. Павлов

Институт динамики геосфер РАН, Москва, 119334, Россия

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

Ключевые слова: межблоковые перемещения, динамические воздействия, расчетная модель

Numerical simulations of interblock displacement accumulation under low-amplitude dynamic action

A.M. Budkov, G.G. Kocharyan and D.V. Pavlov

Institute for Dynamics of Geospheres RAS, Moscow, 119334, Russia

The paper is concerned with numerical realization of a computational model for the effect of seismic vibrations on the stressed contact between rock blocks. The numerical simulations demonstrate that nonlinearity of stress-strain relations and pronounced asymmetry of loading and unloading provide a possibility for accumulation of residual displacements in interblock contact zones under dynamic actions of even very low amplitude compared to the strength of geomaterial. It is found that backward motion in elastoplastic interaction of local contact spots under unloading plays a vital part in the process.

Keywords: interblock displacements, dynamic action, computational model

1. Введение

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

того, представления о строении разломных зон и протекающих там деформационных процессах получаются из данных in situ путем инверсии сейсмологических и других непрямых данных, как правило, в предположении линейной упругости и других упрощений [2].

Любопытно, что на первом месте среди сформулированных в докладе рабочей группы американских сейсмологов Национальному научному фонду США десяти основных задач сейсмологии на ближайшие 1020 лет стоит «большой вызов № 1: Как скользят разломы?» [3]. Более того, не до конца изучены многие детали даже столь, казалось бы, привычного эффекта, как трение между телами. Механика впервые обнаруженного свыше 200 лет назад перехода от статического к динамическому трению до сих пор не вполне понята [4].

© Будков А.М., Кочарян Г.Г., Павлов Д.В., 2010

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

Одной из первых удачных компьютерных моделей для исследования различных, в том числе динамических, процессов деформирования дискретных сред была модель Кандела [5], в которой горный массив представлялся жесткими частицами — дисками, а деформации локализовывались в точках контактов. С тех пор появилось множество компьютерных реализаций дискретных моделей, например коды 3DEC [6], DDA [7], метод клеточных автоматов [8, 9] и т.д. При этом все расчеты для дискретной среды демонстрируют исключительную важность правильного описания закономерностей деформирования области контакта между блоками.

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

Понятно, что соотношения, пригодные для описания деформационных характеристик межблоковых контактов, могут быть получены лишь из эксперимента. Выполненные лабораторные опыты позволили в деталях исследовать закономерности динамического деформирования трещин различных типов [10], оценить влияние скорости де формирования [11] и значимость процессов восстановления прочности [12]. Проведение экспериментов со взрывами разного масштаба дало возможность определить характеристики разломов и трещин в условиях естественного залегания и сопоставить закономерности деформирования натурных объектов с образцами лабораторного масштаба [1].

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

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

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

2. Данные экспериментов

В Институте динамики геосфер РАН в течение ряда лет сейсмическими методами проводились инструментальные исследования характеристик межблоковых границ разных иерархических уровней — от мелких трещин до региональных разломов. Результаты измерений деформационных свойств нарушений сплошности приведены, например, в работах [1, 10]. Там же детально изложена и методика наблюдений. Из полученных данных следует, что отличительной особенностью реологических зависимостей нарушений сплошности массива горных пород является нелинейность связи между действующим напряжением и деформацией разлома или трещины.

На рис. 1 приведены динамические диаграммы «напряжение - деформация» для модельной трещины толщиной около 2 мм, заполненной песком [1]. Можно видеть, что по мере роста деформации сдвиговая жесткость трещины ^ = ёт/ёп постепенно снижается. В этом соотношении т — касательное напряжение; и — относительное перемещение бортов трещины. Заметим, что в области малых деформаций ход различных диаграмм полностью совпадает, что свидетельствует о достоверности полученных результатов. Ниспадающие

Рис. 1. Диаграммы «напряжение - деформация» для различных уровней динамического воздействия на межблоковый контакт, заполненный песком. Пунктирная линия соответствует наклону 15 000 МПа/м

Относительная деформация, микрострейн

Относительная деформация, микрострейн

Рис. 2. Изменение средней жесткости разломной зоны (а) и модельной трещины (б) с ростом амплитуды динамического воздействия

участки диаграмм (ветви разгрузки) имеют один и тот же наклон —15 000 МПа/м, что соответствует величине сдвиговой жесткости контакта при малых деформациях. Подчеркнем, что жесткость динамической разгрузки контакта остается неизменной при разном уровне воздействия на контакт.

В ряде случаев бывает удобно использовать понятие «средней» жесткости, линеаризуя ветвь нагружения:

к (1) п(/тах)

где ^шах — момент времени, в который напряжение достигает максимума.

Как показывают результаты измерений [1, 10, 14] (рис. 2), величина как нормальной, так и сдвиговой средней жесткости снижается с ростом максимальной деформации у тах в соответствии с эмпирическим соотношением

к =-------ко----------------------------------. (2)

1+(у пи*/ у» Г

При этом величина к0 соответствует жесткости разгрузки, а параметр у» изменяется в диапазоне 10-9-10-6 в зависимости от напряженного состояния и типа нарушения сплошности [1, 10]. При сдвиговом деформировании зависимость более сильная (в соотношении (2) показатель степени т — 0.6-0.8), чем при нормальном (т — 0.2-0.3).

Вполне определенной является и зависимость средней жесткости от масштаба нарушения сплошности. Если исключить из рассмотрения мелкие трещины длиной менее 100 м, то данные измерений с коэффициентом корреляции R = 0.98 описываются регрессионной зависимостью

к = 837L_0'41, (3)

где жесткость k измеряется в МПа/м, а длина разлома L — в километрах [15].

Типичная реологическая зависимость «сдвиговое напряжение - перемещение» схематично показана на

рис. 3 сплошной линией. Характерными точками этой диаграммы являются пиковое перемещение пр — точка достижения предельной прочности на сдвиг тр, а также смещение пг, при котором прочность выходит на остаточное значение тг. Эти параметры, как и отношение тр1 тг, изменяются в довольно широких пределах и зависят от таких факторов, как шероховатость контакта, прочность поверхности, свойства материала-заполнителя, Р-Г-условия и т.д. При этом величина пр составляет 10-3-10-2 от характерного размера сдвигаемого блока, а выход на остаточное значение прочности происходит при очень больших сдвиговых смещениях пг — 100 пр.

Характеристики контакта не остаются неизменными в процессе деформирования. При этом может наблюдаться как увеличение, так и снижение жесткости меж-блокового контакта, а следовательно, и скорости накопления межблокового перемещения. Известно, что и сила сопротивления сдвигу, и жесткость контакта зависят от скорости относительного перемещения берегов трещины, времени стационарного контакта и амплитуды перемещения [16].

ир и,

Сдвиговое перемещение

Рис. 3. Типичная реологическая зависимость для сдвигового деформирования контакта между блоками горной породы

Рис. 4. Изменение сдвиговой жесткости трещины к8 при многократном динамическом воздействии; т/тр = 0.69 (а), 0.74 (б); 1-4 — этапы деформирования

Увеличение жесткости обусловлено двумя факторами. Первый, последовательное увеличение жесткости контакта при многократных циклических нагружениях, известен для квазистатического циклического нагружения [14]. При динамическом нагружении трещин этот эффект выражен слабее, однако также довольно существенен. На рис. 4, а показано изменение сдвиговой жесткости трещины при многократном динамическом воздействии. Описание методики проведения экспериментов, в которых исследовалась динамика деформирования межблоковых контактов, и подробное изложение результатов можно найти в работах [10-13]. Как видно из рис. 4, а, величина к8 постепенно возрастает, приближаясь к значению жесткости разгрузки, которое показано на графике пунктирной линией. Здесь же сплошной линией показано соотношение

k = -

k0

1 + 8.9N

-0.63

(4)

в котором к0 = 15 000 МПа/м — жесткость разгрузки.

Если деформация контакта происходит со скоростью ниже некоторого критического значения, то жесткость увеличивается пропорционально логарифму времени [11]. В наших экспериментах критическое значение скорости относительного смещения блоков составляло 0.05-0.1 мкм/мин или 25-50 мм/год, что замечательным образом совпадает с верхней границей характерных значений скорости смещения по сейсмогенным разломам [11]. При слишком больших скоростях деформации упрочнения нарушения сплошности не происходит, а следовательно, маловероятно и возникновение динамической неустойчивости.

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

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

[18].

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

3. Расчетная модель

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

- Взаимодействие динамического импульса с нарушением сплошности в напряженной среде может привести к возникновению остаточного межблокового перемещения из-за нелинейности соотношений «напряжение - деформация» для нарушений сплошности. Необратимые деформации накапливаются при циклическом воздействии.

- Основным макроскопическим критерием возникновения неустойчивости является достижение напряженно-деформированным состоянием контакта спадающей ветви реологической зависимости.

- Последовательное увеличение жесткости контакта при многократном воздействии приводит к «привыканию» среды к уровню динамической нагрузки.

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

Для построения расчетной модели зависимость жесткости от перемещения на восходящей ветви зависимости «напряжение - деформация» зададим аналогично (2):

к = к

s0

, )а

1 — u

up

u<u

p

(5)

В этом случае восходящая ветвь кривой описывается соотношением

1—

uju

du = ks0up I(1 — x)adx =

ks0up

a +1

1 —

\a+1

(6)

Для учета эффекта увеличения жесткости при повторных воздействиях, в выражении для жесткости (5) показатель степени а зависит от количества циклов N.

Аналитическое описание запредельного участка зависимости т(и) существенно зависит от типа контакта и Р-Г-условий. Для определенности, далее будем использовать соотношение

1

Т = IТ p

1 + exp

1—

У/

(7)

При этом если текущее состояние контакта (п0, т0) не соответствует зависимости (6), т.е. кумулятивное перемещение в силу истории процесса деформирования выше значения, определяемого (7), то жесткость нагружения определяется зависимостью (5) с соответствующей заменой п0 на п0 - Ап (рис. 5). Если кумулятивное перемещение превышает значение пр, то жесткость нагружения снижается по сравнению с (5) в соответствии с изменением отношения тр1 /Тр , где Тр — прочность контакта, а тр1 — текущее значение прочности. Значение жесткости разгрузки к80 постоянно в ходе всего процесса деформирования.

Для проведения численных расчетов использовался двумерный вычислительный код [19], разработанный на основе лагранжева численного метода «Тензор» [20].

1р

тр1

То

/

1

-I J

/ All 1 ! і 1 i ' i i

и0 ир

Перемещение Рис. 5. Схема к построению расчетной модели

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

Уравнения, описывающие движение и напряженное состояние твердого деформируемого материала в осесимметричном случае, имеют вид:

■dP+p div v = 0, dt

Vr =■

dr , dt ,

Vz =-

dz , dt ,

P dvr дsrr дsrz 2srr + Szz + дР

p---------—-----------—-----------—---------------------+---------= 0,

dt дг дz r дг

P

srz дР

— + ^~ = g, г дz

(8)

dvz дszz _дІ1_ dt дz дг

dє

Р~ _ sггeгг _ szzezz _ s00e00 _ dt

_2s Є _Р — = 0 р dt

где t — время; г, z, 0 — цилиндрические координаты ^ — ось симметрии); р — плотность; , Vz — компоненты вектора скорости V; g — ускорение свободного падения; Р — давление; sij — девиатор тензора напряжений; Єу — девиатор тензора скоростей деформаций; є — удельная внутренняя энергия; — лагранжева

производная по времени: //dt = д//дt + (V, V)/.

Для плоского случая в декартовой системе координат уравнения (8) можно записать в виде:

dp

dt

+ р div v = 0,

dx dy

v x =---, v y = —,

x dt y dt

P-

dt дx ду дx

(8а)

dv у дs дs дР

Р—y----------------------------yyy-il +-= g

dt ду дx ду

dє

I-------

dt

P — — sxxexx — syyeyy — szzezz —

— 2s e —— — = 0 xy xy p dt '

a

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

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

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

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

4. Результаты расчетов и обсуждение

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

критерия прочности необходимым условием возникновения остаточного перемещения при квазигармоничес-ком динамическом воздействии на контакт будет условие превышения суммой статических т8( и динамических касательных напряжений предела прочности на-

рушения тр:

Т* +т<1 >тр- (9)

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

т = С + ст (10)

где C — сцепление; ф — угол трения на поверхности раздела; стп — эффективное нормальное напряжение. Задача решалась в плоской двумерной постановке в декартовой системе координат. Упругие характеристики материала подвижного блока взяты типичными для скальных пород: плотность р0 = 2.7 г/см3, скорость продольных волн а0 = 5 км/с, скорость поперечных волн Ь0 = 2.55 км/с. Для более отчетливого проявления эффектов накопления перемещений было задано довольно малое значение коэффициента трения на поверхности скольжения tg ф = 0.1, а сцеплением для простоты было решено пренебречь: C = 0. Эти упрощения не имеют принципиального значения.

На первом этапе расчета к верхней грани блока прикладывалось статическое напряжение стп и решалась задача на установление. Во всех расчетных вариантах полагалось стп = 0.1 МПа. Для ускорения процесса сходимости при проведении этого этапа расчета были существенно (в 2-10 раз) увеличены коэффициенты искусственных вязкостей и каждые 100 шагов обнулялись скорости во всех узлах расчетной сетки.

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

р =@р0 ^0 при г^ (11)

1 @р0 при * > *0-

Амплитуда Р0 варьировалась в окрестности фрикционной прочности контакта т = ст^ф = 0.01 МПа. Время нарастания нагрузки г0 принято равным 10 мс. Проведенные расчеты показали, что поведение рассматриваемой модели в случае статической нагрузки вполне соответствует классическим представлениям. Если сдвигающая сила превышает силу трения, то происходят срыв блока и его равноускоренное движение. Если же сдвигающая сила меньше силы трения, наблюдаются

4 6

Время, мс Время, мс

Рис. 6. Эпюры смещения и скорости движения правого нижнего угла блока

слабые колебательные движения блока и его деформация под действием сдвигающей силы. Движения блока как целого при этом не происходит.

Значительно больший интерес представляют численные эксперименты, в которых при допредельном смещающем усилии (Р0 = 0.009 МПа) в блоке дополнительно возбуждались упругие колебания. Динамическая нагрузка моделировалась воздействием на верхнюю грань блока треугольного импульса давления с амплитудой вдвое меньше нормального статического давления Рк = 0.05 МПа и длительностью гк = 1 мс. Приложенное импульсное воздействие возбуждает в блоке колебания, которые приводят к изменению нормального напряжения на границе скольжения и, соответственно, силы трения. В фазах разгрузки сила трения становится меньше сдвигающего усилия, что приводит к возникновению периодического проскальзывания блока. Этот процесс иллюстрирует рис. 6, на котором показаны эпюры смещения и скорости в правой нижней угловой точке блока, где дискретный характер движения проявляется наиболее отчетливо. По мере затухания колебаний фазы проскальзывания становятся все короче. В приведенном варианте расчета продолжительность относительного движения составила около 15мс.

Использование критерия (9) приводит, казалось бы, к ясной картине механики накопления межблоковых перемещений. Подобная схема использована, например, в модели «акустической флюидизации» Дж. Мелоша [21]. Эксперимент между тем демонстрирует несостоятельность такой модели. Во-первых, накопление перемещений происходит при величинах , значительно меньших, чем это следует из соотношения (9) [13]. Во-

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

В основной серии расчетов мы использовали изложенную выше (раздел 3) модель, в которой возможность накопления межблоковых деформаций на напряженном контакте при малых амплитудах динамического воздействия связана с нелинейностью реологических характеристик нарушений сплошности. В численных экспериментах рассматривалась система из двух блоков по конфигурации близкая к проведенным ранее лабораторным опытам [1, 14]. Подвижный блок А высотой h = = 3.2см и длиной L = 8 см расположен на стержне В длиной 1 м (рис. 7). Межблоковый контакт характеризуется следующими основными константами модели Тр = 0.04 МПа, тг = 0.02 МПа, к80 = 800 МПа/м, пр = = 0.01 см. В начальный момент времени на левой боковой грани подвижного блока в соответствии с соотношением (11) задавалось сдвигающее усилие в виде граничного давления РЦ. Амплитуда нагрузки Р0 = = 0.06 МПа, время нарастания г0 =1 мс. Поскольку отношение длины блока к высоте Цк = 2.5, то для равновесия необходимо, чтобы на контакте блоков «мобилизованная» величина напряжения сопротивления сдвигу составляла значение т» = 0.024 МПа, что значительно меньше заданной в расчетах величины максимального сцепления тр, однако больше остаточной прочности контакта тг.

Рис. 7. Постановка задачи (вторая серия расчетов)

Время, мс

Рис. 8. Эпюры горизонтальных составляющих скорости относительного

После затухания переходных процессов (в момент времени ? = 20 мс) к левому торцу стержня прикладывается имитирующая удар импульсная нагрузка Рк (г) с амплитудой 0.5 МПа и длительностью 0.1 мс.

На рис. 8 приведены полученные в расчете временные эпюры скорости и перемещения середины нижней грани подвижного блока относительно основания. На расчетных эпюрах отчетливо видны колебания, сопровождающие процесс установления равновесия после приложения сдвигающей нагрузки РЦ. Эти колебания обусловлены наличием упругой составляющей сдвиговой деформации межблокового контакта, которая, как будет показано ниже, играет важную роль в исследуемом процессе. На рис. 9 приведена зависимость «напряжение - межблоковое перемещение», рассчитанная для всего процесса деформирования. В соответствии с используемой моделью в результате асимметрии цикла «нагрузка - разгрузка» при наличии сдвигающего усилия образуется незамкнутая гистерезисная петля, т.е. возникает остаточное межблоковое перемещение. При многократном динамическом воздействии, которое в расчете смоделировано упругими колебаниями, возбуждаемыми ударом по торцу стержня в момент времени 1 = 20 мс, точка, соответствующая напряженно-деформированному состоянию контакта в текущий момент времени, перемещается вправо вдоль оси абсцисс (рис. 9). На рис. 8 этому процессу соответствует постепенное нарастание величины перемещения.

После того как напряженно-деформированное состояние контакта достигнет спадающей ветви реологической зависимости, сопротивление сдвигу начинает

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

Проведенные численные эксперименты продемонстрировали, что нелинейная реология приводит к воз-

Время, мс

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

можности кумуляции остаточных перемещений на меж-блоковом контакте. Однако следует подчеркнуть, что в основе отклика блочной системы на малые деформации лежит не только, или, вернее, не столько кулоновское трение между поверхностями блоков, сколько упругопластические деформации межблокового контакта. Ку-лоновское трение, в его классическом виде, подразумевает, что сила сопротивления сдвигу всегда направлена в сторону противоположную перемещению. При этом если движения нет, то сила трения всегда направлена против сдвигающей силы и равна ей по величине. В численном эксперименте мы можем провести расчет, наложив ограничения на контактное взаимодействие:

» [-Т^п V при V ф 0,

ТГ = I п (12)

?- т8 при V = 0,

где т* — «скорректированное» сопротивление сдвигу;

— величина сопротивления сдвигу, рассчитанная в соответствии с нелинейной моделью (раздел 3); т8 — локальное сдвигающее напряжение, обусловленное внешней динамической или статической нагрузкой.

0.00 -|----1---1----1----1----1---1----1----1----1---г-

0 100 200 300 400 500

Относительное перемещение, мкм

Рис. 9. Траектория изменения напряженно-деформированного состояния межблокового контакта. Сплошная линия — результаты расчета, пунктирная — реологическая зависимость, определяемая соотношениями (6), (7)

10

20 30

Время, мс

40

50

Время, мс

Рис. 10. Эпюры горизонтальных составляющих скорости относительного движения и относительного перемещения середины основания блока

Подобные расчеты были проведены в той же постановке, что и описанные выше (рис. 7). Результаты такого расчета приведены на рис. 10.

Как и следовало ожидать, после приложения сдвигового усилия относительное движение уже к моменту времени 1 = 2 мс останавливается после достижения величины относительного перемещения и = 68 мкм, которая соответствует восходящей ветви реологической зависимости (п). После возбуждения в системе упру-

гих колебаний (1 > 20 мс) начинается медленное прерывистое движение блока, аналогичное показанному на рис. 6, которое имеет место лишь в случае выполнения условия (9). Полученный результат еще раз наглядно показывает, что в предложенной модели эффект накопления деформации связан не только с нелинейным характером зависимости т(и), но и с возможностью возникновения возвратного движения («отдачи») при разгрузке контакта.

5. Заключение

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

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

Разработанный математический аппарат и программное обеспечение могут быть эффективно использова-

ны для численного моделирования процессов динамического деформирования блочного горного массива.

Работа выполнена при финансовой поддержке РФФИ (гранты №№ 09-05-12023-офи-м и 09-05-00968-а).

Литература

1. Кочарян Г.Г., Спивак A.A. Динамика деформирования блочных массивов горных пород. - М.: ИКЦ «Академкнига», 2003. - 423 с.

2. Li Y.-G., Vidale J.E., Aki K., Xu F. Depth-dependent structure of the Landers fault zone from trapped waves generated by aftershocks // J. Geophys. Res. - 2000. - V. 105. - No. B3. - P. 6237-6254.

3. Seismological Grand Challenges in Understanding Earth’s Dynamic Systems: Report to the National Science Foundation, IRIS Consortium / Ed. by T. Lay. - 2009. - 76 p.

4. Rubinstein S.M., Cohen G., Fineberg J. Dynamics of precursors to frictional sliding // Phys. Rev. Lett. - 2007. - V. 98. - No. 2. - P. 226103 (4 pages).

5. Cundall P.A. A Computer Model for Simulating Progressive Large Scale Movements in Blocky Rock Systems // Proc. Int. Symp. on Rock Fracture, Paper II-8. - Nancy: ISRM, 1971. - P. 129-136.

6. Jing L., Hudson J.A. Numerical methods in rock mechanics // Int. J. Rock Mech. Min. Sci. - 2002. - V. 39. - No. 4. - P. 409-427.

7. Shi G.H. Three-Dimensional Discontinuous Deformation Analysis // Proc. 4th Int. Conf. on Discontinuous Deformation Analysis (ICADD-4), Glasgow, June 6-8, 2001. - Glasgow: University of Glasgow Scotland, UK, 2001. - P. 1-22.

8. Псахье С.Г., Остермайер Г.П., Дмитриев A.M., Шилько Е.В., Смолин AM., Коростелев С.Ю. Метод подвижных клеточных автоматов как новое направление дискретной вычислительной механики. 1. Теоретическое описание // Физ. мезомех. - 2000. - Т. 3. - №2.-С. 5-15.

9. Aстафуров С.В., Шилько Е.В., Псахье С.Г. Влияние стесненных условий на характер деформирования и разрушения блочных сред при сдвиговом нагружении // Физ. мезомех. - 2009. - Т. 12. - N° 6.-

С. 23-32.

10. Костюченко В.Н., Кочарян Г.Г., Павлов Д.В. Деформационные характеристики межблоковых промежутков различного масштаба // Физ. мезомех. - 2002. - Т. 5. - № 5. - С. 23-42.

11. Кочарян Г.Г., Кулюкин A.A., Mарков В.К., Mарков Д.В., Пер-ник ЛЖ. О критической скорости деформации разломных зон // ДАН. - 2008. - Т. 418. - № 3. - С. 383-386.

12. Кочарян Г.Г., Павлов Д.В. Нарушение и залечивание зон локализации деформаций в массиве горных пород // Физ. мезомех. -2007.- Т. 10. - № 1. - С. 5-18.

13. Кочарян Г.Г., Кулюкин A.A., Павлов Д.В. Роль нелинейных эффектов в механике накопления малых возмущений // Физ. мезомех. -2006. - Т. 9. - № 1. - С. 5-14.

14. Кочарян Г.Г., Костюченко В.Н., Павлов Д.В. Инициирование деформационных процессов в земной коре слабыми возмущениями // Физ. мезомех. - 2004. - Т. 7. - № 1. - С. 5-22.

15. Кочарян Г.Г. Физический смысл отклонения некоторых параметров сейсмического процесса от закона подобия // ДАН. - 2009. -Т. 429. - № 6. - С. 821-824.

16. Dieterich J.H. Time-dependent friction and the mechanics of stick-slip // Pure Appl. Geophys. - 1978. - V. 116. - No. 4-S. - P. 790-806.

17. Scholz C.H. The Mechanics of Earthquakes and Faulting. - New York: Cambridge Univ. Press, 1990. - 439 p.

18. Кочарян Г.Г., Павлов Д.В. Макроскопический критерий возникновения фрикционной неустойчивости на контакте между блоками горной породы // Физические поля и динамика взаимодействую-

щих геосфер: Сб. научных трудов ИДГ РАН. - М.: ГЕОС, 2007. -С. 93-104.

19. Aрхипов В.Н., Борисов ВЛ., Будков A.M. и др. Механическое действие ядерного взрыва. - М.: Физматлит, 2002. - С. 50-58.

20. Mайнчен Дж., Сак Е. Метод расчета «ТЕНЗОР» // Вычислительные методы в гидродинамике. - М.: Мир, 1967. - С. 185-211.

21. Melosh H.J. Acoustic fluidization: A new geologic process? // J. Geophys. Res. - 1979. - V. 84. - No. B13. - P. 7513-7520.

22. Кочарян Г.Г., Кулюкин A.A. Модель накопления малых динамических возмущений в блочном горном массиве // Физические поля и динамика взаимодействующих геосфер: Сб. научных трудов ИДГ РАН. - М.: ГЕОС, 2007. - C. 83-92.

Поступила в редакцию 17.02.2010 г.

Сведения об авторах

Будков Александр Михайлович, д.т.н., снс ИДГ РАН, jack77@mail.ru Кочарян Геворг Грантович, д.ф.-м.н., зав. лаб. ИДГ РАН, gevorgk@idg.chph.ras.ru Павлов Дмитрий Вячеславович, к.ф.-м.н., внс ИДГ РАН, pavlov@idg.chph.ras.ru