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

УДК 539.3 + 550.34

Экспериментальное исследование различных режимов скольжения блоков по границе раздела. Часть 3. Численное моделирование

A.M. Будков, Г.Г. Кочарян

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

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

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

Experimental study of different modes of block sliding along interface.

Part 3. Numerical modeling

A.M. Budkov and G.G. Kocharyan

Institute of Geosphere Dynamics RAS, Moscow, 119334, Russia

This paper completes a series of studies on the laws of block sliding along interface. It has been shown that in order to model the whole range of crustal block movements, the empirical rate-and-state friction law must be supplemented with a term that accounts for the appearance of additional shear resistance associated with the dynamic viscosity of the contact between blocks. With this term, the experimentally observed slow sliding events can be modeled with good accuracy. A generalization of results of the entire series of studies published in several issues of the journal suggests that both the dynamic and quasi-static modes of sliding along faults are components of a single deformation process. The parameter that governs the formation and evolution of a fault sliding mode is the ratio between the effective values of the fault zone rigidity and rock mass rigidity. Their variation determines the occurrence of a particular sliding mode.

Keywords: fault slip modes, rate and state friction law, stick-slip motion, apparent viscosity, numerical modeling

1. Введение

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

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

© Будков A.M., Кочарян Г.Г., 2016

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

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

Во второй публикации [3] были представлены результаты натурных экспериментов, в которых исследовались различные режимы гравитационного соскальзывания блока по естественной поверхности разлома. Было показано, что тип деформационных событий, которые могут быть реализованы, определяется как мезо-структурой берегов трещины, так и вещественным составом ее заполнителя. При деформировании трещины, заполненной неоднородным материалом, миграция областей концентрации напряжений («контактных пятен» или «asperities») в процессе сдвига приводит не только к изменению параметров деформирования, но и к трансформации самого режима вследствие изменения реологических свойств локальных участков межблокового контакта.

Для моделирования и объяснения различных явлений, связанных с трением по границам блоков земной коры (зарождение землетрясений, постсейсмическое скольжение, форшоки, афтершоки), широко используется подход, основанный на экспериментальных исследованиях Дж. Дитериха [4, 5], согласно которому коэффициент трения зависит от скорости скольжения и эволюционирует в процессе перемещения. Эти представления легли в основу «rate and state» (R&S) модели (закон трения Дитериха), в соответствии с которой сила сопротивления сдвигу на единицу площади контакта между блоками представляется в виде

Ts = = an

+ a ln

(i&i ^

u*

V

+ b ln

( uJ6}

D

(i)

где стп — средние эффективные нормальные напряжения; ¡л—текущая величина коэффициента трения; л0 — коэффициент трения, соответствующий стабильному скольжению со скоростью и*; и — скорость перемещения; а, Ь, Dc — эмпирические константы; 6 — переменная состояния, которая определяется из кинетического уравнения

I X | 6

В состоянии стабильного скольжения с некоторой скоростью x = us переменная состояния 6 принимает значение 6 = Dc/us. Это значение используется в качестве начального условия при решении уравнения (2).

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

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

2. Численное моделирование

При численном моделировании для описания фрикционного взаимодействия блока и подложки был использован закон R&S в виде (1), (2). Уравнение, описывающее движение подвижного блока модели, имеет вид

mx = k(ust - x) - Sts, (3)

где m — масса; S—площадь блока; k—жесткость пружины; us — скорость протяжки пружины; Fs = S Ts — сила сопротивления сдвигу, задаваемая в соответствии с уравнением (1).

Для численного решения дифференциального уравнения (3) совместно с уравнениями (1) и (2) использовался метод Рунге-Кутты 4-го порядка аппроксимации.

При воспроизведении в расчетах результатов лабораторных экспериментов параметры R&S-модели задавались так, чтобы обеспечить режим скоростного разупрочнения блока, т.е. А = b - a > 0. При этом значения констант a, b, А и Dc варьировались.

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

d^ dT

dx dx

> K.

(4)

1 = 1 —

D

(2)

При малых положительных значениях величины А изменение параметра Дс может привести к появлению режима условно стабильного скольжения, при этом критические значения параметров в точке бифуркации задается соотношением [6]

А, = ДК. (5)

Соотношение (5) позволяет оценить область изменения параметров, при которых реализуются различные режимы скольжения (стабильный, условно стабильный и прерывистый). Примеры зависимостей скорости движения блока от времени приведены на рис. 1. На этом и следующих рисунках скорость движения блока и нормирована на скорость протяжки и8, а время — на постоянную времени т = Л/ т/К, где т — масса блока; К — жесткость пружины.

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

При неизменной величине параметра А увеличение абсолютных значений а приводит к росту скорости движения блока. При достаточно больших значениях а ы/ы5 ~л[а (рис. 2, а), а длительность фазы скольжения остается практически неизменной. Как известно [4], параметр а определяет величину скачка коэффициента трения d| ~ а в момент резкого изменения скорости скольжения. Параметр Ь определяет степень «разупрочнения» контакта в фазе релаксации. Увеличение диапазона изменения коэффициента трения приводит к увеличению максимальной скорости движения блока и, соответственно, выделяющейся при динамическом срыве упругой энергии. С ростом параметра А = Ь - а амплитуда динамических срывов также растет (рис. 2, б). По мере удаления от критической величины Ас зависимость и(А) стремится к линейной. Длительность пика скорости при этом незначительно сокращается.

Параметр Dc влияет на динамику процесса «разупрочнения» и характеризуется как критическое смещение, при котором коэффициент трения переходит к новому значению. При малых значениях Dc влияние этого параметра на амплитуду динамического срыва довольно слабое (изменение Дс на 3 порядка вызывает изменение величины максимальной скорости в фазе динамического срыва менее чем в 3 раза). Увеличение Дс приводит к росту времени возникновения первого динамического срыва и, соответственно, к увеличению его амплитуды. По мере того как величина параметра Дс приближается к критическому значению, определяемому соотношением (5), начинается переход от режима прерывистого скольжения к режиму условно стабильного скольжения, и амплитуда пиков скорости скольжения снижается тем быстрее, чем ближе Дс к критической величине. В приведенном на рис. 2, в примере рас-

1200 ^т

и/и8' 10010 1 0.1 0.01

0

40

—I—I—Г

1790 1830

3590 3630

5390 5430

^т

Рис. 1. Зависимость скорости блока от времени. Параметры расчетной модели: ы5 = 0.1 мм/с, т = 10 кг, Fn = 100 Н, К = = 1000 Н/м, |0 = 0.4, Дс = 10 мкм, а = 0.001. Режим стабильного скольжения, Ь = 0.0011 (а); режим условно стабильного скольжения, Ь = 0.00119 (б); режим прерывистого скольжения, Ь = 0.009 (в)

чета критическое значение параметра Дс = 53.2 мкм. При этом смена режима скольжения, проявляющаяся в снижении амплитуды пиков скорости смещения блока, начинается при значениях Дс более 10 мкм, т.е. ~20 % от критической величины.

3. Сопоставление с экспериментальными данными

Путем подбора значений констант а, Ь и Дс в расчетах обеспечивалось наилучшее соответствие экспериментальным зависимостям.

и/щ 5000 ; 3000 ■ 2000 ■

1000^

500 : 300 -200 -

100 ^

и/щ :

1000:

100 :

10

0.001 0.01 А

и/щ -¡—I гв]

1000 - !

• • !

500 ■ |

зоо ° ;

200 - |

100 -I—.........—.........—.........—■ ■ ■ .= ■ ■ ■ ■

0.01 0.1 1 10 Дс, мкм

Рис. 2. Зависимость максимальной скорости смещения от параметров закона R&S. Параметры расчетной модели: и8 = = 0.02 мм/с, т = 130 кг, Fn = 1300 Н, К = 16 700 Н/м, ¡0 = = 0.4, скорость движения блока и нормирована на скорость протяжки и8. а — зависимость от параметра а: А = 7 • 104 (1), 1.5 • 103 (2), 2 • 103 (3); б — зависимость от параметра А, вертикальный пунктир соответствует критической величине А с = 2.6 • 10-4, а = 0.01 (1), 0.001 (2), 0.0001 (3), и8 = 2 • 10-5 м; в — зависимость от параметра Дс, а = 0.01, Ь = 0.0107, вертикальный пунктир соответствует критической величине Дс = = 53.2 мкм.

На рис. 3 сопоставлены расчетная и экспериментальная зависимости скорости скольжения блока от времени для контакта, заполненного кварцевым песком. Для удобства здесь и далее время отсчитывается от момента максимума скорости.

В расчете использовались следующие параметры R&S-модели: а = 0.0002, Ь = 0.00109, Дс = 10 мкм. Характеристики «слайдера-модели соответствовали эксперименту. Можно видеть, что в случае выраженного прерывистого скольжения удается достаточно точно вос-

Рис. 3. Зависимость скорости скольжения от времени для контакта, заполненного кварцевым песком. Параметры расчета: К = 17 КН/м, л = 0.61, Дс = 10 мкм, а = 2 • 10-4, Ь = 1.09 • 10-3, П = 0, А с = 1.3 • 10-4. Сплошная линия — расчет, пунктир — эксперимент

произвести экспериментальную эпюру. Сопоставление результатов расчета и эксперимента на диаграмме усилие-перемещение также демонстрирует хорошее соответствие. Расчеты, проведенные для множества экспериментов, показали, что эпизоды прерывистого скольжения («лабораторные» землетрясения) удовлетворительно моделируются при помощи канонического закона R&S (1).

Попытки численно смоделировать режимы медленного скольжения (лабораторные «медленные землетрясения» и «эпизоды медленного скольжения») натолкнулась на непреодолимые сложности. Это хорошо видно из примера, приведенного на рис. 4, где показаны результаты эксперимента и численного расчета эпизода скольжения блока по контакту, заполненному сухой глиной при жесткости нагружающего элемента 17 кН/мм. В эксперименте время скольжения в опыте с сухой глиной возрастает по сравнению с контактом из кварцевого песка примерно на порядок, при соответствующем сни-

Рис. 4. Зависимость скорости скольжения от времени для контакта, заполненного сухой глиной: 1 — эксперимент, 2 — расчет, Па = 0; 3 — расчет, = 1000 Па • с. Параметры модели: жесткость пружины К = 17 кН/м, А с = 6.5 • 104, л = 0.7, Дс = = 50 мкм, а = 1.2 • 10-4, Ь = 1.35 • 10-3

Время ^ - ¿тах, с Время ^ - ¿тах, с

Рис. 5. Зависимость скорости скольжения от времени: а — контакт, заполненный обводненной глиной, 1 — эксперимент, 2 — расчет, параметры модели: К = 10 кН/м, Ас = 6.5 • 104, | = 0.56, Дс = 90 мкм, а = 1.0 • 104, Ь = 3.5 • 10-3, п = 2.76 • 104 Па • с; б — контакт, заполненный смесью кварцевого песка с тальком, 1 — эксперимент, 2 — расчет, параметры модели: К = 16.56 кН/м, А с = 2.3 • 10-4, | = 0.61, Дс = 18 мкм, а = 1.0 • 10-5, Ь = 8.5 • 10-4, Па = 3.94 • 104 Па • с

жении и максимальной, и средней скорости скольжения (кривая 1 на рис. 4).

Сигнал такого типа не удается воспроизвести в расчетах с использованием канонической R&S-модели (линия 2 на рис. 4). Варьируя параметры, можно подобрать амплитуду сигнала, но наполненность эпюры, т.е. величина смещения, в основном определяется массой блока и жесткостью пружины (параметры, строго определенные для конкретного эксперимента). В расчете возможно лишь повлиять на «размывание» переднего фронта пика, изменяя параметр Дс. Расчетный пик динамического срыва при заданной амплитуде неизменно остается существенно более узким, чем полученный в эксперименте.

Результаты численного моделирования «медленных» движений удалось приблизить к эксперименту путем введения в каноническое уравнение R&S (1) члена, учитывающего возникновение дополнительного сопротивления сдвигу, связанного с динамической вязкостью контакта:

|0 + а 1п

+ Ь 1п

м*0 Д

Па

" d '

(6)

где Па — коэффициент динамической вязкости межблокового контакта; d—площадь и толщина контактной зоны.

Результаты расчета с учетом динамической вязкости на рис. 4 представлены линией 3. Можно видеть, что ввод в уравнение движения члена, связанного с динамической вязкостью контакта, позволяет в расчете с удовлетворительной точностью воспроизвести и «медленные» движения.

Использование уравнения (6) позволило путем подбора эффективной вязкости удовлетворительно воспроизвести характер движения во всех экспериментах. На рис. 5 приведены примеры эпюр скоростей смещения

блока в двух опытах — с контактом, заполненным увлажненной глиной, и с контактом, заполненным кварцевым песком с добавлением 25 % талька.

В этих случаях движение имеет длительную фазу (~30-40 с) плавного увеличения скорости смещения, а затем фазу торможения примерно той же продолжительности. Значение максимальной скорости смещения снижается до 10-100 мкм/с. Для наилучшего совпадения результатов расчета и эксперимента коэффициент «эффективной вязкости» п был увеличен более чем на порядок по сравнению с предыдущим расчетом.

4. Обсуждение и выводы

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

Зависимость величины коэффициента эффективной вязкости полученного путем достижения наилуч-

Рис. 6. Зависимость величины максимальной скорости смещения блока от отношения жесткостей контакта и пружины. Жесткость пружины К= 10.8 (1), 16.9 (2), 30.9 Н/мм (3)

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

П = ^и,^^ (7)

описывающей данные с коэффициентом корреляции Д = 0.95.

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

Рис. 7. Зависимость величины эффективной вязкости от максимальной скорости смещения блока

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

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

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

Работа выполнена при поддержке РНФ (проект № 14-17-00719).

Литература

1. Кочарян Г.Г., Новиков B.A. Экспериментальное исследование различных режимов скольжения блоков по границе раздела. Часть 1. Лабораторные эксперименты // Физ. мезомех. - 2015. - Т. 18. -№4. - C. 94-104.

2. Kocharyan G.G., Kishkina S.B., Novikov V.A., Ostapchuk A.A. Slow slip events: parameters, conditions of occurrence, and future research prospects // Geodynam. Tectonophys. - 2014. - V. 5(4). - P. 863-891.

3. Кочарян Г.Г., Остапчук A.A., Павлов Д.В., Ружич B.B., Батух-тин И.В., Виноградов E.A., Камай A.M., Марков В.К. Экспериментальное исследование различных режимов скольжения блоков по границе раздела. Часть 2. Полевые эксперименты и феномено-

логическая модель явления // Физ. мезомех. - 2015. - Т. 18. - № 6. -С. 75-85.

4. Dieterich J. Modeling of rock friction: 1. Experimental results and constitutive equations // J. Geophys. Res. B. - 1979. - V. 84. - No. 5. -P. 2161-2168.

5. Popov VL., Grzemba B., Starcevic J., Popov M. Rate and state dependent friction law and prediction of earthquakes: What can we learn from laboratory models? // Tectonophysics. - 2012. - V 532. - P. 291-300.

6. Scholz C.H. Earthquakes and friction laws // Nature. - 1998. - V. 391.-

P. 39-42.

7. Белинский И.В., Михалюк A.B., Христофоров Б.Д. Вязкость горных

пород при деформационных процессах // Физика Земли. - 1975. -№ 8. - C. 80-84.

8. Ци Ч., Ван М., Цянь Ц., Чень Ц. Структурная иерархия и механичес-

кие свойства горных пород. Часть I. Структурная иерархия и вязкость // Физ. мезомех. - 2006. - Т. 9. - № 6. - C. 29-39.

9. Христофоров Б.Д. Исследование реологических характеристик твердых тел в широком диапазоне времен деформирования // Физ. мезомех. - 2010. - Т. 13. - № 3. - С. 111-115.

10. Ide S, Shelly D.R., Beroza G.C. Mechanism of deep low frequency earthquakes: Further evidence that deep non-volcanic tremor is generated by shear slip on the plate interface // Geophys. Res. Lett. - 2007. -V. 34. - P. L03308. - doi 10.1029/2006GL028890.

11. Peng Z., Gomberg J. An integrated perspective of the continuum between earthquakes and slow-slip phenomena // Nature Geosci. -2010. - V. 3. - P. 599-607. - doi 10.1038/ngeo940.

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

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

Будков Александр Михайлович, д.т.н., снс ИДГ РАН, [email protected] Кочарян Геворг Грантович, д.ф.-м.н., зам. дир. ИДГ РАН, [email protected]