CC BY
66
Моделирование поведения высокопористых геоматериалов при формировании полос локализованного уплотнения
Ю.П. Стефанов12, М. Тьерселен2
1 Институт физики прочности и материаловедения СО РАН, Томск, 634021, Россия 2 Московский научно-исследовательский центр Шлюмберже, Москва, 109004, Россия
В работе предложены модификации моделей Друккера-Прагера-Николаевского и Рудницкого и построены алгоритмы, необходимые для численного моделирования процессов локализованного уплотнения высокопористых геоматериалов. Данные модели могут быть использованы для описания процессов после уточнения параметров и согласования результатов численного и натурного экспериментов. С целью отладки и проверки их работоспособности смоделировано поведение слоя уплотняющегося материала в условиях сжатия при различных граничных условиях и рассмотрено влияние некоторых параметров на результаты расчетов. Показано, что предложенные модели позволяют моделировать формирование и распространение полос локализованного сжатия в высокопористых средах.
Modeling of the behavior of highly porous geomaterials at compaction band formation
Yu.P. Stefanov1,2 and M. Thiercelin2
1 Institute of Strength Physics and Materials Science SB RAS, Tomsk, 634021, Russia
2 Schlumberger Moscow Research Center, Moscow, 109004, Russia
In the paper we put forward modifications of the Drucker-Prager-Nikolaevskii and Rudnicki models. The algorithms necessary for the numerical modeling of compaction band formation in highly porous geomaterials are constructed. These models can be used to describe processes after parameters are adjusted and the results of numerical and full-scale experiments are correlated. In order for the model efficiency to be adjusted and verified, we model the behavior of a compacted material layer under compression in different boundary conditions and study the influence of certain parameters on the calculation results. It is shown that the considered models allow one to model the formation and propagation of compaction bands in highly porous media.
1. Введение
В геологических средах выделяют три основных типа неустойчивого поведения, которые определяются особенностями структуры, механическими свойствами материала и условиями нагружения. Такая классификация типов неустойчивого поведения и условий их возникновения приведена в работе [1]. Это — локализация сдвига, уплотнения и расширения, а также смешанные процессы, включающие явления, характерные для смежных типов поведения. Локализация сдвига — вероятно, наиболее часто наблюдаемая форма неустойчивости, которая обычно возникает в плотных консолидированных средах, а также в высокопористых мате-
риалах при умеренном давлении обжатия или среднего напряжения. В процессе развития полос локализованного сдвига часто имеет место разрыхление среды, сопровождающееся объемными изменениями, или дилатан-сией. Полосы разуплотнения, расширения соответствуют отрывному механизму разрушения.
Формирование полос локализованного уплотнения остается наименее изученной формой проявления неустойчивости в геологических средах. Данное явление наблюдается в высокопористых средах (обычно более 15 %) при больших значениях среднего давления (более 100 МПа) [1-7]. Причем, если процесс формирования полос сдвига нередко сопровождается некоторым рас-
© Стефанов Ю.П., Тьерселен М., 2007
ширением, ростом пористости, то при образовании полос уплотнения происходит значительное ее снижение. Образование полос (зон) уплотнения сопровождается сильным измельчением частиц, переходом в катаклас-тический режим деформирования и резким снижением проницаемости.
Экспериментальное исследование поведения высокопористых образцов [2-7] показало, что неоднородная пористая структура материала не выдерживает высоких давлений, происходит дробление частиц, в результате чего их обломки забиваются в поры, снижается сопротивляемость сдвиговой нагрузке. Формируются узкие полосы уплотнения, часто ориентированные перпендикулярно оси наибольшего сжатия. Разрушение частиц возникает в результате того, что напряженно-деформированное состояние на масштабе частиц неоднородно даже при внешнем гидростатическом сжатии. Изучение структуры материала в областях локализации при запредельном давлении показало существенное снижение пористости.
Однако интерес представляет не поведение отдельных частиц, тем более, рассмотрение поведения всех частиц невозможно, а поведение макроскопической среды или образца. Поэтому описание процесса деформирования осуществляется на уровне усредненных макроскопических моделей. В работах [3-6] построены предельные поверхности для ряда материалов. Анализ экспериментальных данных показал, что поверхность, описывающая прочность высокопористых геологических сред, в осях интенсивности напряжений и среднего напряжения имеет замкнутый вид шапки, что говорит о снижении прочности по мере увеличения величины среднего напряжения после определенного значения. В связи с этим основные модели, описывающие данное поведение, основаны на поверхностях такого замкнутого вида cap-модели. В настоящее время для описания поведения высокопористых материалов с учетом формирования зон уплотнения предложен ряд моделей [8-14]. Учитывая характерный вид предельных поверхностей, построенных на обобщении экспериментальных наблюдений, уравнения всех моделей содержат член с квадратичной зависимостью от среднего напряжения (давления).
Следуя теории пластического течения, неупругое (пластическое) поведение материала определяется в соответствии с заданными предельной поверхностью и законом течения:
f (о, е, %) = 0,
dep = ,
где f— поверхность текучести; % — параметр упрочнения; g — пластический потенциал; X — множитель, определяемый в ходе процесса; ер — компоненты пластической (неупругой) деформации; а. — компоненты
тензора напряжений. Последнее выражение можно записать в несколько ином виде [8]:
dep- = dA, ij
Sn 1
— Gx — GCTS j т T 3 J
V J
где sij = Ор -1/3 акк8р- — компоненты девиатора напряжений; — компоненты тензора напряжений; а = = -акк/3 — первый инвариант напряжений; т = = (sijSij| 2)1/2 — второй инвариант девиатора напряжений; Gx =дg/дт; Ga =дg/да. Интенсивность пластической деформации сдвига dYр = 2(ёе| ёер / 2)12 запишется в виде dYр = GтdЛ, а объемная часть пластической деформации ёер = -ёеР : ёер = GadA, где Яр — компоненты тензора деформаций; ер =Е р « 8 р/3.
Рудницкий предложил использовать для описания процесса деформирования с учетом уплотнения предельную поверхность эллиптической формы [8-11] с полуосями, длины которых есть функции накопленных неупругих объемной и сдвиговой деформации. Таким образом, по мере развития пластической деформации происходит изменение формы эллипса в соответствии с характером упрочнения. Эллиптическую форму поверхности предлагают Ди Маджио и Сандлер [13], а также де Борст и Гроен [14]. Модель Керрола [12] предполагает подвижную параболическую поверхность, Сван и Сео [15] предложили вариант комбинированной поверхности, состоящей из трех частей: поверхности Кулона-Мора, круговой части в области высокого давления и ограничения в области растяжения.
Таким образом, все предложенные уравнения предельных поверхностей для описания формирования полос уплотнения включают квадратичную зависимость от давления. К сожалению, в приведенных работах не уделяется внимания пластическому потенциалу, без которого невозможно осуществить вычисления. В частности, большинство авторов допускают применение ассоциированного закона, т.е. в качестве поверхности потенциала использовать предельную поверхность. Однако это существенно сужает возможности вычисления, т.к. означает жесткую привязку коэффициента дилатансии к предельной поверхности, что не позволит удовлетворить экспериментальным данным, а также приведет к дополнительным сложностям и ограничениям диапазона нагрузки. Безусловно, использование неассоциированного закона потребует построения необходимой функции, для чего необходимо осуществить дополнительные экспериментальные измерения.
При использовании уравнения поверхности в форме Рудницкого (рис. 1) [8-11]:
г = (°-г/ +4-1.
2 7.2
а Ь
Для удобства запишем уравнение в виде:
/ = (а-Ро)2 + R^2 - а2, где р0 — центр, а, Ь — полуоси эллипса; R = а2/Ь2. Подобную форму уравнения примем и для пластического потенциала
g = (ст-Ро)2 + Ят2,
пластическая деформация и напряжения будут определяться из выражений:
ёуР = 2^, ёеР = 2Х(а- ро), х*
X =
а =
2^ +1’ а* + 2ХКро
1 + 2ХК
Параметр X может быть найден из решения уравнения:
А2 ( т* А2
+ Я
в - Ро
1 + 2ХК
2ХБц +1
-а2 = 0.
Для модели Керрола [12]
I =
(Р-Ро)
+ -
1
а~ Ь
так же перепишем уравнение:
I = (ст-Ро)2 + Ят2 - а2 Здесь ро, а, Ь — параметры параболы; Я = а2/Ъ.
При пластическом потенциале в виде:
g = (ст-Ро)2 + Эх, мы имеем следующие выражения для вычисления деформаций и напряжений:
dYр =А£,
dєp = -2Х(а- Ро),
% = %* - АЭц,
_ а* + 2ХКро
~ 1 + 2ХК ’
а параметр X определяется из уравнения:
-Ро
V
1 + 2ХК
+ Я(х* - АЭД ■
= о.
Успехом в изучении особенностей формирования полос локализованного уплотнения является формулировка условий их образования и возможной ориентации к осям нагружения [16-18]. Например, в случае осесимметричного нагружения условия формирования полосы уплотнения приобретают вид: "П + Р < -л/3", где "Л = = (д// да)/(д// Эх) — коэффициент наклона касательной к предельной поверхности и Р = -ёер /ёу р — коэффициент дилатансии.
Несмотря на активные исследования последних лет в области изучения особенностей поведения высокопо-
ристых геологических материалов и формирования полос уплотнения в литературе практически не встречаются результаты численного моделирования таких процессов. В работе [19] выполнены численные расчеты по моделированию поведения образца с кольцевым вырезом в условиях одноосного сжатия, однако несмотря на использование соответствующей модели формирование полосы локализованного уплотнения не наблюдалось. В работах [20, 21] рассматривается распространение участка уплотнения и сопоставляется напряженное состояние образца, содержащего отрезок уплотнения, с состоянием в образце с трещиной отрыва.
2. Процедура расчета
Для моделирования процесса деформации будем численно решать систему уравнений механики, включающую уравнения движения и неразрывности, а в некоторых случаях, при определенных формах уравнения состояния, и уравнение энергии. Эта система замыкается определяющими соотношениями, которые конкретизируют поведение среды, устанавливая связь между напряжениями, деформациями и/или скоростями, определяют ее поведение, а также геометрические соотношения для вычисления тензора деформации. Краткое изложение используемого подхода для описания процессов деформации можно найти в работах [22-24]. Рассмотренные ниже задачи решались в двумерной постановке для условий плоской деформации. Решение системы уравнений осуществлялось с использованием явной конечно-разностной схемы из [25].
Таким образом, описание процесса деформации будем осуществлять путем численного решения системы уравнений механики сплошной среды:
^ у,р _ ’
Ро*о =Р^.
Здесь ро, р — начальное и текущее значение плотности материала; У0, V — начальное и текущее значение элементарного объема материала.
Будем считать, что полные деформации £у состоят из упругой е®- и пластической ер- частей. Аналогичное разложение примем и для скорости деформации:
= е ® + ер у у у.
Связь между напряжениями и деформациями для упругого поведения опишем гипоупругим законом:
у =-а5 у + sij,
°'1] = 2^(еу -1Е 1к5у I,
Dt
Dt
' - Яу Яїк тук Яук ,
а = - К
V V'
Здесь Еу — компоненты тензора деформаций Коши: 1,
=~ (Ч у + у),
иг- — компоненты вектора перемещения; Юу — компоненты тензора скоростей вращения:
= '2(иг-,у ~ иу, ^
точка сверху означает производную по времени; К и ^ — модули сжатия и сдвига соответственно.
На каждом последовательном шаге по времени приращения пластической деформации будут пропорциональны разнице между напряжениями, рассчитанными по упругому закону, и напряжениями, соответствующими предельной поверхности (текучести). Первый этап расчета после определения координат точек расчетной сетки и расчета полных приращений деформации состоит в предварительном расчете напряжений по упругому закону, обозначая при этом предварительно вычисленные напряжения звездочкой:
(уП+1)* = уП + (ДуП+1)*,
(ап+1)* = ап + (Дап+1)*.
После этого проверяется, находится ли данная точка в пространстве напряжений внутри предельной поверхности или нет, т.е. проверяется условие перехода в пластическое состояние. Проверка условия осуществляется подстановкой значения напряжений в уравнение, определяющее предельную поверхность. Если /* = = /(о*-) < о, то состояние материала в данной ячейке расчетной сетки в текущий интервал времени находится внутри предельной поверхности, т.е. данный элемент среды находится в упругом состоянии и можно продолжить расчет на последующий слой времени, опустив символ *. Другими словами, рассчитанное напряженное состояние соответствует истинному и можно переходить к дальнейшему расчету на следующем слое по времени.
В случае, если /* = /(Оу) > о, т.е. точка в пространстве напряжений оказалась вне поверхности, то условие
пластичности считается выполненым: элемент среды перешел в пластическое состояние и часть деформации имеет неупругий характер. Необходимо вычислить эту часть пластической деформации, так чтобы точка в пространстве напряжений оказалась на поверхности, определяющей область упругого состояния и ограничивающей возможное напряженное состояние среды.
Использование ассоциированного закона означает, что из данного предварительно вычисленного напряженного состояния точка 2* должна быть перенесена на предельную поверхность, переместившись к ней по нормали (рис. 1, 2). В случае использования неассоциированного закона течения, когда уравнения предельной поверхности и пластического потенциала не совпадают, точка в пространстве напряжений перемещается по нормали к поверхности пластического потенциала до пересечения с предельной поверхностью. Таким образом осуществляется корректировка напряжений, чтобы удовлетворить условиям, определенным моделью поведения.
Предельная поверхность и поверхность пластического потенциала могут изменяться от шага к шагу деформирования, поверхности эволюционируют от одного временного шага к другому в соответствии с заданными законами.
Для осуществления расчетов, моделирующих поведение среды в ограниченной расчетной области, нам необходимо получить выражения, которые позволят вычислить напряжения и деформации в ходе процесса.
3. Модели поведения материала
Рассмотрим две модели: первая представляет собой модифицированную модель Друккера-Прагера-Нико-лаевского [26-29] с неассоциированным законом течения и повреждениями, ограничивающими состояние в области сжатия, вторая — модифицированную модель Рудницкого [8-11].
Рис. 1. Эллиптическая предельная поверхность модели Рудницкого: а и Ь — полуоси, р0 — центр эллипса
Рис. 2. Предельная поверхность в модифицированной модели Друкке-ра-Прагера-Николаевского. р0Х и Р02 — давление начала и полной деградации материала; фи у — углы, соответствующие внутреннему трению и дилатансии; X* — точка, соответствующая предварительному расчету напряжения по упругому закону
3.1. Модифицированная модель Друккера-Прагера-Николаевского с неассоциированным законом течения и ограничением состояния в области сжатия
В качестве уравнения предельной поверхности возьмем уравнение, соответствующее по форме уравнению Друккера-Прагера [26]:
/ = т-аа- с.
Допустим, что параметр а, имеющий физический смысл коэффициента внутреннего трения, зависит от величины среднего напряжения:
а = а(ер, а) = ао(ур)(1 - П(а, ер)),
/ _ л2
где В(а) =
Ро1
для ро1 < о < ро2, аналогич-
ч ро2 ро1 у ным образом запишем выражение для с:
с = с(ер, а) = Со (ур )(1 - В(а, ер)).
Таким образом, с, а, определяющие предельную поверхность, не остаются постоянными, а меняются как в ходе деформирования, так и с ростом давления. Используемый метод рассматривает малые приращения напряжений и деформаций и допускает выполнение расчетов при изменении их значений на каждом временном интервале вычислений по мере изменения состояния материала. Параметры, определяющие поверхность течения и направление вектора пластической деформации, приобретают смысл функций от давления, накопленной сдвиговой и/или объемной пластической деформации, значения которых постоянны лишь в течение одного временного шага. При достижении давлением в точке значения, когда начинаются разрушение зерен и уплотнение среды, происходит рост функции, которой можно придать смысл деградации. Причем, при разгрузке (снижении давления) и при последующем нагружении мы автоматически будем следовать поверхности с учетом предшествующих процессов деградации среды.
Подставляя в уравнение поверхности текучести выражения для коэффициентов с учетом деградации, получим:
1-
в- рог ро2 _ ро1
ап + х-
1-
ро1
ро2 _ ро1
Таким образом, мы перешли от линейной предельной поверхности к поверхности более сложного типа. Получив падающий характер предельной поверхности, мы сохранили смысл коэффициента а как угла внутреннего трения. Особенностью такого описания является то, что вид поверхности будет отличаться от линейного лишь при активном нагружении с ростом давления. Еще одной особенностью данной модели является то, что не возникает проблемы сопряжения поверхностей различного типа, характерной для комбинированных моделей.
Параметры ро1 и ро2 также сделаем изменяющимися в зависимости от накопленной объемной пластической деформации, что обеспечит расширение области упругого состояния и упрочнение среды:
ро1 = ро1 + rЛРо2,
1 + И
1-
Р
У
где Ро, Р — начальный и текущий коэффициенты дила-тансии.
Вид полученной предельной поверхности показан на рис. 1.
Для вычисления неупругой деформации необходимо определить направление вектора ее приращения в пространстве напряжений, в котором строится предельная поверхность. С этой целью воспользуемся линейной формой пластического потенциала, приняв неассоциированный закон течения, как это было предложено Николаевским [27-29]. Получим необходимые для вычислений соотношения аналогично работам [22, 23]:
g = -Ро + т.
Тогда приращения пластической деформации и напряжения будут определяться из выражений:
ёу Р =л^ = х,
Эх
аер =Х^ = -Лр,
Эо
т = т* -Ац,
а = а* + Арк,
Х= / */(ц + аКр).
Имея значения второго инварианта девиатора напряжений, вычислим компоненты напряжений как
*у = р
Теперь необходимо ввести ограничение на значения параметров. Очевидно, что для получения однозначного решения необходимо, чтобы существовала точка пересечения нормали к поверхности потенциала и предельной поверхности. Значит, угол между ними не должен приближаться к 90°. Кроме того, величина объемной пластической деформации уплотнения ограничена величиной пористости материала и скорость дилатансии должна снижаться по мере приближения к предельному значению. Примем для скорости дилатансии следующее выражение:
А”
Р = Ро
Необходимо также согласовать упрочнение поверхности и расчет объемной пластической деформации. Очевидно, что скорость упрочнения, т.е. расширения
предельной поверхности, не может превышать определенную величину, определяемую процессом деформации:
Лро2
К
- + Дер < Де.
Максимальная скорость расширения поверхности:
Дро! =^£ = Де,
К К
что выполняется при нулевом значении коэффициента дилатансии, т.е. когда мы фактически размыкаем предельную поверхность вдоль оси давления. Таким образом, когда скорость дилатансии станет равной нулю, поверхность перестанет быть замкнутой, уплотнение прекратится и возможен неограниченный рост давления при дальнейшем уже упругом сжатии вещества, находящегося в «жидком», катакластическом состоянии, пористость которого больше не меняется.
3.2. Модифицированная модель Рудницкого
Рассмотрим вопрос использования предельной поверхности эллиптического типа, предложенной Рудницким [8-11] (рис. 2). Уравнение поверхности имеет вид:
22 а2 Ь2
-1 = о.
Перепишем его для удобства в виде:
R(a-ро)2 + х2 -Ь2 = о, где R = (Ь/а)2.
Как и в предыдущей модели, для выполнения формальных операций, которые позволят получить необходимые для вычислений соотношения, необходимо задать уравнение пластического потенциала. Воспользуемся линейной формой пластического потенциала, принимая во внимание, что использование данного подхода не препятствует изменению значений параметров в ходе процесса деформации, так что все функции строятся для текущего момента времени и в каждой точке среды. Поэтому угол наклона прямой может быть функцией накопленной деформации и/или величины давления. Причем, удобство линейного уравнения состоит в наличии единственного параметра, определяющего наклон поверхности, а соответственно направление вектора пластической деформации. Это существенно упрощает как вычисления, так и подбор параметров для удовлетворения экспериментальным данным:
g = -Ра + т,
ёт р =1^ = 1,
Эх
ёЕр = Х^ = -а.р,
Эо
т = т* -Ац,
а = а* + ЛРК.
Значение параметра X для текущего состояния элемента среды определяется из уравнения:
А2 А + ХВ + С = о,
А = ц2 + ^|ЗК )2,
В = 2(R(a* - ро)рК -х»,
С = /* = /(ст*) = R(a* - ро )2 +х*2 - Ь2,
'.±4в2
^1,2 _ '
■4 АС
2А
Из полученных значений выбирается наименьший положительный корень. Подставляя его в выражения для расчета значений напряжений и деформаций, мы получим напряжения в соответствии с эллиптической предельной поверхностью и значения деформации точки среды в текущий момент времени.
Отдельного рассмотрения требует выбор значений коэффициента дилатансии и скорости изменения (расширения и/или смещения) предельной поверхности в ходе деформирования, которые определяет скорость уплотнения. Очевидно, что скорость неупругой деформации не может оказаться выше, чем скорость полной деформации. Если для сдвиговой части это условие выполняется автоматически и она физически не ограничена, то для объемной деформации необходима формулировка специальных условий. В режиме разрыхления коэффициент дилатансии изменяется в ходе деформирования и стремится к нулю при приближении к состоянию полной деградации среды, когда дальнейшее расширение невозможно или теряет смысл рассмотрение процесса без введения и учета трещин.
В случае уплотнения описание процесса становится более сложным, поскольку снижение модуля отрицательного коэффициента дилатансии до нуля означает фактическое прекращение процесса уплотнения. Причем дальнейший рост давления не ограничен. Среда переходит в катакластическое состояние, в котором невозможно лишь растяжение без потери сплошности, но в условиях сжатия она ведет себя как жидкость. Напряженное состояние не может быть вне предельной поверхности, а использование замкнутой поверхности приводит к определенным условия и для сжатия. Для преодоления этого необходимо изменение поверхности по мере сжатия. Мы говорим лишь о сжатии. Допустим, что сдвиговая прочность уже исчерпана и напряженное состояние лежит на гидростатической оси. Фактически, после потери сдвиговой прочности, когда состояние материала оказалось на гидростатической оси, наша предельная поверхность уже превратилась в отрезок, лежащий на оси гидростатического сжатия, а его длина ограничивает величину давления. Таким образом, скорость изменения длины этого отрезка определяет, сопровождается ли процесс деформирования неупругой
объемной деформацией, т.е. уплотнением. Теперь скорость эволюции поверхности становится однозначно связана со скоростью дилатансии, так что нулевое значение коэффициента дилатансии означает расширение поверхности со скоростью деформирования и, наоборот, бесконечное значение коэффициента дилатансии означает, что материал не сопротивляется сжатию, что очевидно некорректно. Поэтому, аналогично тому, как это сделано для предыдущей модели, для коэффициента дилатансии запишем:
р=р.
о
а для параметров, определяющих предельную поверхность:
а = ао + Да,
с = Со + гДа,
Да = -КДе(1 -Р/Ро)п.
Таким образом мы согласовали скорость объемного упрочнения и скорость уплотнения. В начальный момент, когда ер = о, мы имеем Р = Ро и Да = о. По мере уплотнения и уменьшения пористости коэффициент дилатансии снижается, приращение полуоси предельной поверхности увеличивается и при ер ^ е”ах полу-
чаем р ^ о и Да = Да. модуль будет:
Отсюда текущий объемный
к = К
1 -
( Є Р -ЕР
°тах °
4. Результаты расчетов
При изучении поведения геологических материалов, построении предельных поверхностей и зависимостей напряженного состояния от деформации возникает сложная проблема, связанная с тем, что процессы протекают неоднородно. Модель поведения среды строится на основании экспериментальных данных о поведении образцов конечных размеров, при этом предполагается, что поведение каждой точки образца соответствует поведению образца в целом.
Обратимся к проблеме формирования полос уплотнения. В случае осесимметричного нагружения происходит не просто формирование отдельных полос локализации уплотнения, но и изменение формы образца. В зонах локализации материал выдавливается из образца в стороны (например, [2, 7]), часто формируются несколько участков интенсивного развития неупругой деформации и даже в пределах одной полосы локализации степень деформации может быть различной. Напряженное состояние неоднородно не только в пределах всего образца, но и в плоскости локализации. При построении кривых в координатах (Р = а, Q =
= а1 -а°, где а1 — осевая нагрузка; оо —среднее обжимающее давление) в каждой точке образца данная
Рис. 3. Схема нагружения образца с уплотняющимся слоем
величина может различной. Поэтому модель и подобранные параметры не гарантируют, что поведение элемента среды в зоне локализации следует построенным соотношениям, т.к. путь нагружения в каждой точке может быть различным. Однозначную оценку работоспособности модели и ее пригодности для описания поведения не единичного элемента среды, а физического объекта — образца — могут дать лишь результаты численного моделирования.
Чтобы решить вопрос о возможности использования описанных в разделе 3 моделей, пока не касаясь более общей проблемы полного описания поведения образца, смоделируем поведение тонкого слоя уплотняющегося материала, помещенного между блоками высокопрочной среды с теми же упругими характеристиками (рис. 3). В таких условиях напряженно-деформированное состояние в слое компактируемого материала будет однородным, что дает возможность проверки и отладки моделей, рассматривая напряженно-деформированное состояние материла с ростом давления.
Рис. 4. Пути нагружения в тонком слое компактируемого материала в условиях одноосного деформирования при различной величине начального обжимающего давления без учета упрочнения, рассчитанные с использованием модифицированной модели Друккера-Пра-гера-Николаевского
Рис. 5. Пути нагружения в компактируемом слое при различной величине начального давления в условиях одноосного деформирования и соответствующий вид предельной поверхности (а); изменение параметров модели от величины давления в ходе деформирования (б)
Рис. 6. Изменение напряженного состояния в компактируемом слое в ходе деформирования (е — осевая деформация слоя) (а); изменение неупругих сдвиговой и объемной частей деформации с ростом давления (б)
Таблица 1
Параметры модели для расчетов на рис. 5, 6
с, ГПа Рої, ГПа Р02, ГПа а0 Р0 г h п £ Р тах т
0.0025 0.075 0.2 0.5 -0.75 0 1.0 1 0.1 2
На рис. 4-9 показаны результаты расчетов поведения компактируемого слоя материала с использованием модифицированной модели Друккера-Прагера-Нико-лаевского. В соответствии с принятой моделью, параметры которой зависят от величины давления и накопленной неупругой объемной деформации, напряженное состояние (рис. 4, 6) следует предельной поверхности, показанной на рис. 1. Причем в ходе нагружения после достижения предельной поверхности может быть учтено упрочнение среды при перемещении точки предельного давленияр02, принимая во внимание уменьшение пористости (рис. 5). Можно предположить, что в рамках данной модели изменение параметра р02 обеспечит комплексный учет возможного упрочнения среды, т.к. напряженное состояние следует предельным поверхнос-
тям, которые имеют линейный вид с положительным коэффициентом а, зависящим от величины давления.
Хорошо видно, что предлагаемый вид функций (рис. 5, б) изменения параметров обеспечивает переход в неупругое состояние и начало уплотнения вещества в заданной форме предельного состояния, которая соответствует уплотнению материала. Скорость уплотнения зависит от пути нагружения, а значит и от напряженного состояния, при котором начинается неупругое деформирование, а также от предельной величины уплотнения — уменьшения пористости. Параметры модели для результатов, представленных на рис. 5, 6, приведены в табл. 1.
На рис. 6 показаны изменения напряженного состояния в ходе деформации, а также величина объемной и сдвиговой деформации с ростом давления. По мере уплотнения материала меняется наклон кривой зависимости напряжений от деформаций: чем интенсивнее развивается процесс уплотнения, тем более пологий характер кривых. На стадии, когда скорость уплотнения максимальна, хорошо виден горизонтальный участок
Рис. 7. Пути нагружения в различных точках компактируемого слоя при отсутствии кинематических ограничений на боковых гранях образца (а); параметры напряженного состояния в тех же точках (б): 1 — точка в центральной части; 2 — вблизи поверхности образца. Q = ау -а°
Таблица 2
Параметры модели для расчетов на рис. 7
с, ГПа Р01, ГПа Р02, ГПа а0 Р0 г к п ер тах т
0.0025 0.075 0.2 0.3 -0.75 0 1.0 2.0 0.1 0.5
на кривой зависимости осевой компоненты напряжения и давления от деформации.
При возможности поперечной деформации образца, т.е. при отсутствии кинематических ограничений на боковых гранях образца, напряженное состояние в плоскости локализации в ходе развития неупругой деформации становится неоднородным. Это хорошо видно на рис. 7, где показаны путь нагружения и характеристики напряженного состояния для двух точек в уплотняющемся слое. При построении кривых зависимости напряжений от деформации (рис. 7, б) использовались параметры, как это принято при построении экспериментальных зависимостей: по оси абсцисс отложена осевая деформация, а по оси ординат — разница между действующим напряжением и боковой нагрузкой, имеющей постоянное значение. На рисунке хорошо видно различие путей нагружения этих точек среды. В центральной части слоя наблюдается быстрое снижение сдвиговой прочности и больший рост давления, тогда как вблизи поверхности образца после заметного снижения прочности наблюдается ее увеличение. Параметры модели, соответствующие рис. 7, приведены в табл. 2.
Важное значение для формирования полос локализованного уплотнения имеют условия нагружения, которые можно разделить на две составляющие: вид граничных условий и начальная величина обжимающего давления. В соответствии с условиями формирования полос локализованного уплотнения [16-18], чем выше величина начального давления, тем больше предпосылок локализованного характера уплотнения, т.к. в процессе дальнейшего нагружения достижение предельной поверхности происходит в области с большим углом оги-
бающей. Аналогичное влияние имеют граничные условия и вид нагружения. Таким образом, при одноосном деформировании создаются наиболее благоприятные условия формирования полос локализованного уплотнения, поскольку угол наклона пути нагружения в упругой области в этом случае составляет
02 ,
1 + V
тогда как при осесимметричном нагружении 01 = V?. Для процесса формирования полос локализации уплотнения также существенное значение имеет скорость упрочнения среды, на что влияют величина пористости и диапазон ее возможного изменения.
На рис. 8, 9 приведено несколько примеров моделирования формирования полос локализованного уплотнения. В однородном материале наблюдается как образование рассеянных участков локализации, так и формирование периодических полос. Следует отметить, что с ростом нагрузки по мере упрочнения и исчерпания возможности уплотнения полосы локализации будут расширяться и вовлекать в данный процесс остальной материал.
При наличии участка ослабленного материала, в котором инициируется процесс уплотнения, с ростом нагрузки возможны прямолинейное продвижение зоны локализации и формирование полосы. Нельзя забывать, что развитие неупругой деформации, и тем более при локализованном ее характере, является динамическим процессом, сопровождающимся динамическим перераспределением напряжений и излучением упругих волн. Поэтому на развитие данного процесса могут оказать влияние не только квазистационарное напряженное состояние, но и динамические эффекты, как вызванные внешним воздействием, так и иницирован-ные собственным процессом и отражениями от различных границ, приходом волн из одной вершины в другую, аналогично явлениям, характерным для динамического распространения трещин.
I
-1.12630Е—3 -1.12807Е—3 -1.12984Е-3 -1.13161Е-3 -1.13338Е—3 -1.13515Е-3 -1.13692Е—3 -1.13869Е-3 -1.14046Е-3
-1.14223Е-3
-1.14400Е-3
-1.39290Е—2 —2.01764 Е+0 -4.02134 Е+0 -6.02505Е+0 -8.02876Е+0 -1.00325Е+1 -1.20362Е+1 -1.40399Е+1 -1.60436Е+1 -1.80473Е+1 -2.00510Е+1
Ш—
-6.22425Е—1
-9.33383Е—1
-6.39175Е+0
-8.94825Е+0
Рис. 8. Распределение объемной пластической деформации -£Р при формировании полос локализованного уплотнения в последовательные моменты времени при различной величине начального обжимающего давления: а0 = 150 (а); 300 МПа (б)
При очень высоком, близком к предельному, начальном давлении возмущение напряженного состояния может привести к зарождению множества рассеянных участков уплотнения (рис. 9). Ориентация и расположение таких участков будут определяться не только основным напряженным состоянием, но и дополнительным воздействием-возмущением, которое может иметь различную природу, в том числе и связанным с распространением полосы локализованного уплотнения. Такие зоны начального уплотнения могут не диагностироваться экспериментально, т.к. интенсивность неупругой деформации в них имеет существенно меньшую величину, чем в основной зоне локализации, а возможность их развития зависит от дальнейшего роста нагрузки.
и
0.00000Е+0 рп---------
Ш -
-1.10572Е-1 ■ -------
—2.21144Е—1 —3.31716 Е—1 .
—4.42288Е—1 —5.52861 Е—1 -6.63433Е—1 -7.7/005Е 1 -8.84577Е—1 -9.95149Е—1 —1.10572Е+0 "
—1.23991 Е+0 —1.37410Е+0
Аналогичные приведенным выше расчеты были выполнены с использованием модели с эллиптической предельной поверхностью, описанной в разделе 3. На рис. 10, 11 построены кривые изменения состояния с ростом нагрузки при различной величине предварительного сжатия. При отсутствии упрочнения видно, что состояние материала следует неподвижной предельной поверхности эллиптической формы (рис. 10). При учете упрочнения среды начальной предельной поверхности соответствует лишь начало неупругого деформирования. В ходе дальнейшей деформации скорость упрочнения зависит от скорости уплотнения, т.е. коэффициента дилатансии, а также управляющего параметра п, определяющего вид этой зависимости. Влияние параметра п на путь нагружения и, как следствие, на скорость уплотнения показано на рис. 11, 12. Более высокому показателю степени п соответствуют быстрое снижение прочности и короткий пологий участок кривой
Рис. 9. Объемная неупругая деформация -£р (а) и повреждения (б) при распространении полосы локализованного уплотнения при высоком обжимающем давлении
Рис. 10. Пути нагружения в тонком слое компактируемого материала в условиях одноосного деформирования при различной величине начального обжимающего давления (а0 = 0.0, 0.01, 0.05, 0.1 ГПа) при неподвижной эллиптической предельной поверхности
Рис. 11. Пути нагружения в компактируемом слое при различной величине начального давления в условиях одноосного деформирования и соответствующий вид начальной предельной поверхности (п = 1) (а); влияние параметра п на путь нагружения (б)
Рис. 12. Изменения напряженного состояния в компактируемом слое деформации с ростом давления при различных значениях п (б)
Таблица 3
Параметры модели для расчетов на рис. 11, 12
а, ГПа Ь, ГПа Р0, ГПа Р0 ер °шах п т
0.175 0.045 0.085 -0.85 0.2 1.0; 1.5; 2.0 1
ходе деформирования (а) и неупругих сдвиговой и объемной частей
антитрещиной. Данное сопоставление кажется оправданным, т.к. впереди вершины такого отрезка возникает концентрация напряжений, что способствует его продвижению и уплотнению материала вблизи вершины. К этому следует добавить, что специфика изменения свойств материала в зоне уплотнения также создает дополнительные предпосылки такого продвижения. В
зависимости напряжений от деформации (рис. 12, а) и более высокой скоростью уплотнения (рис. 12, б). Параметры модели, соответствующие результатам на рис. 11, 12, приведены в табл. 3.
Использование данной модели, также как и рассмотренной выше модификации модели Друккера-Прагера-Николаевского, позволяет описать формирование полос локализованного уплотнения (рис. 13). На характер и периодичность этих полос большое влияние оказывают начальный коэффициент дилатансии (рис. 13), а также величина начального давления и условия нагружения, поскольку от этих условий зависит путь нагружения и соответственно точка на предельной поверхности, где начинается процесс уплотнения.
На рис. 14 показано напряженное состояние в образце с распространяющейся полосой локализованного уплотнения. Влияние отрезка локализованного уплотнения нередко сопоставляют с трещиной отрыва, или
-1.36300Е—1 -1.03394Е+0 -1.93158Е+0 —2.82922Е+0 —3.72686Е+0 —4.62450Е+0 —5.52214Е+0 -6.41978Е+0 —7.31742Е+0 -8.21506Е+0 —9.11270Е+0
—1.18230Е—2 —8.87791Е—1 —1.76376Е+0 —2.63973Е+0 —3.51569Е+0 I -4.39166Е+0 I —5.26763Е+0 I -6.14360Е+0 I —7.01956Е+0
I —7.89553Е+0 I -8.77150Е+0
Рис. 13. Объемная пластическая деформация -£р при множественном формировании полос локализованного уплотнения: р0 = - 06
( а ) ; - 0.2 (б)
I
п 0.00000Е+0 |б ■ ■ 2.40810Е—1 |в —1.10520 Е—1 0 11111111111 ■
| —1.10970Е+0 1 1 —1.24431Е—1 Щ
[ —2.21940Е+0 1 ■ 2.2 6616Е—1 -1.31387 Е—1 —1.39733Е—1 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! I
—3.32910Е+0 1 1 2.124 21Е—1 -1.48080Е-1 -1.56427Е-1 III!!!!!!!!!!!!!!!! 1
1 2.06880Е—1 Ц 2.01341 | ■ -1.64774Е-1 11
—5.54850Е+0 Ц || 1.9579! ! Щ -1.73120 Е—1 | 11111 11 I
-6.65820Е+0 8 1 1.90258Е 1 1 1.84717Е-1 —1.89814Е—1 1 |||||||||:: -
—7.76790Е+0 Ц 1.79176Е-1 —1.98160 Е—1 !!!!!!!!!!:||||||!!!111111!!!1 III
1 -8.87760Е+0 1 1 1.73635Е-1 ::: 1.68094Е-1 ШШ1ШШШШШШ Н —2.06507 Е—1 —2.14854Е—1 II
В —9.98730Е+0 ::: 1.62553Е-1 1.57012Е—1 1 -2.27832Е-1
1.53250Е-1 —2.40810 Е—1
-1.95420Е—1
-2.00207Е-1 -2.04072Е-1 -2.07938Е-1 -2.11803Е—1 -2.15668Е—1 -2.19533Е—1 -2.23399Е—1 -2.27264Е-1 —2.31129Е—1 -2.34994Е—1 —2.38860Е—1
Рис.14. Распределения объемной неупругой деформации -Ер (%) (а) и напряжений (ГПа) а (б); ах (в); ау (г) при распространении полосы локализованного уплотнения
т, ГПа 0.06
0.04 -
0.02
0.00
0.0
2
<' ■— 1
/ / » / / ч / / *• / / : \ Г 5 '
0.1
0.2
0.3
а, ГПа
Рис. 15. Пути нагружения в различных точках компактируемого слоя, а также для условий одноосной деформации (пунктирная линия) (а), параметры напряженного состояния в тех же точках (б): 1 — точка в центральной части; 2 — вблизи поверхности образца. Q = ау - а0
Таблица 4
Параметры модели для расчетов на рис. 15
а, ГПа Ь, ГПа Po, ГПа Ро Ер °тах п т
0.175 0.045 0.085 -0.85 0.2 2.0 0.5
зоне уплотнения происходит снижение сдвиговой прочности, постепенно материал переходит в катакласти-ческое состояние, напряженное состояние приближается к гидростатическому. Из условий равновесия следует перераспределение напряжений, приводящее к росту давления в данной зоне, тем самым оказывая дополнительное давление впереди вершины, в еще не охваченной процессом области.
На рис. 15 показаны кривые, иллюстрирующие изменение напряженного состояния с ростом деформации и давления аналогично рис. 7, построенные для двух точек в плоскости уплотняющегося материала в условиях одноосного приращения нагрузки, т.е. без ограничения горизонтальной деформации при постоянном значении боковой нагрузки. Пути нагружения в этих точках совпадают лишь на упругой стадии деформирования и по мере развития неупругой деформации все сильнее расходятся. Параметры модели, соответствую-
< -
У
—7.57813Е+0
-8.06250Е+0 -8.54688Е+0 —9.03125Е+0 —9.51563Е+0 -1.00000Е+1
Рис. 16. Распределение давления (ГПа) (а) и неупругой объемной деформации -£р (б) в образце с уплотняющимся слоем
щие расчетам на рис. 15, приведены в табл. 4. Рисунок 16 иллюстрирует процесс деформации при локализованном уплотнении, где хорошо видны различия напряженного состояния и деформации в пределах уплотняющегося слоя.
Исходя из графиков на рис. 15, б и рис. 7, б, становится понятно, что при построении экспериментальных зависимостей рассмотрение усредненной величины Q далеко не всегда обеспечивает необходимую информацию для получения необходимых функций и предельных поверхностей.
5. Заключение
В исследовании формирования полос локализованного уплотнения к настоящему времени для большого набора материалов экспериментально определен диапазон напряжений, при котором подобные полосы формируются. Построены предельные поверхности. Предложены сар-модели для описания уплотнения высокопористых материалов. Определены условия, при выполнении которых происходит формирование полос локализованного сдвига или полос компактирования, представляющие собой соотношения для коэффициентов упрочнения, наклона огибающей предельной поверхности и коэффициента дилатансии, а также напряженного состояния.
В данной работе рассмотрены модификации моделей Друккера-Прагера-Николаевского и Рудницкого, которые, по мнению авторов, могут быть использованы для описания процесса локализованного уплотнения геоматериалов. Построены необходимые для проведения численных расчетов алгоритмы, и выполнен анализ влияния параметров на развитие процесса. На ряде примеров показано, что предложенные модели позволяют описывать процессы локализованного уплотнения, однако для описания поведения конкретных материалов необходима дальнейшая работа по анализу экспериментальных зависимостей, а также выбору функций и уточнению параметров моделей.
Рассмотрение предложенных вариантов моделей показало, что в некоторых случаях обе модели приводят к идентичным результатам, но, в целом, более устойчивым алгоритмом обладает модель с эллиптической поверхностью, особенно при описании поведения материала на этапе уплотнения с высокой скоростью. Это обусловлено ограничением на соотношение коэффициентов внутреннего трения и дилатансии в модифицированной модели Друккера-Прагера-Николаевского.
Модель с эллиптической поверхностью обеспечивает более простой выбор значений параметров для описания расширения поверхности. В то же время, ряд особенностей модели линейного типа обеспечивает более понятную интерпретацию в привычных терминах механики геоматериалов. Кроме того, применение модели такого типа кажется более оправданным для учета
деградации материала при описании повторного нагружения. Вероятно, будет полезным построение комбинации этих моделей, используя за основу модель с эллиптической поверхностью при расчете напряжений и деформаций, но рассматривая изменения состояния и предельной поверхности с привлечением представлений модели Друккера-Прагера-Николаевского.
Предложенные модели позволяют описывать формирование и распространение полос локализованного сжатия высокопористых средах. Показано, что напряженное состояние, создаваемое перед вершиной полосы локализованного уплотнения, способствует ее продвижению, что соответствует характеру антитрещины отрыва, некоторым особенностям, связанным с ростом гидростатической части напряжений в полосе. Численное моделирование позволило показать, что напряженно-деформированное состояние при развитии деформации в полосе локализованного сжатия неоднородно, соответственно не совпадают в различных точках среды и пути нагружения, что создает определенные сложности в анализе экспериментальных данных.
Авторы выражают благодарность В.Н. Николаевскому и И.А. Гарагашу за полезные дискуссии по проблемам неупругой деформации геоматериалов и численного описания процессов локализованного уплотнения.
Литература
1. Schultz R.A., Siddharthan R. A general framework for the occurrence and faulting of deformation bands in porous granular rocks // Tectono-physics. - 2005. - V. 411. - P. 1-18.
2. Olsson W.A. Theoretical and experimental investigation of compaction bands in porous rock // J. Geophys. Res. - 1999. - V. 104. -P. 7219-7228.
3. Cuss R.J., Rutter E.H., Holloway R.F. The application of critical state soil mechanics to the mechanical behaviour of porous sandstones // Int. J. Rock Mech. Mining Sci. - 2003. - V. 40. - P. 847-862.
4. Schutjens P.M.T.M., Hanssen T.H., Hettema M.H.H., Merour J., deBreeP., Coremans J.W.A., Helliesen G. Compaction-induced porosity permeability reduction in sandstone reservoirs: Data and model for elasticity-dominated deformation // SPE Reservoir Evaluation and Engineering. - 2004. - V. 7. - No. 3. - P. 202-216.
5. Bieda A.El., Sulema J., Martineau F. Microstructure of shear zones in Fontainebleau sandstone // Int. J. Rock Mech. Mining Sci. - 2002. -V. 39. - P. 917-932.
6. Karner S.L., Chester F.M., Kronenberg A.K., Chester J.S. Subcritical compaction and yielding of granular quartz sand // Tectonophysics. -
2003. - V. 377. - No. 3-4. - P. 357-381.
7. Fortin J., Stanchits S., Dresen G., Gueguen Y Acoustic emission and velocities associated with the formation of compaction bands in sandstone // J. Geophys. Res. - 2006. - V. 111. - B10203. - doi:10.1029/ 2005JB003854.
8. Holcomb D.J., Rudnicki J.W Inelastic constitutive properties and shear
localization in Tennessee marble // Int. J. Numer. Anal. Meth. Geo-mech. - 2001. - V. 25. - P. 109-129.
9. Grueschow E., Rudnicki J.W. Elliptic yield cap constitutive modeling for high porosity sandstone // Int. J. Solids Struct. - 2005. - V. 42. - P. 4574-4587.
10. Rudnicki J.W. Compaction Bands in Porous Rock in Bifurcations and Instabilities in Geomechanics // Proceedings of the Int. Workshop on Bifurcation and Instability, St. John’s College, Minnesota, June 3-5, 2002 / Ed. by J. Labuz, A. Drescher. - Swets & Zeitlinger (formerly Balkema) Publishers, 2002. - P. 29-39.
11. Rudnicki J.W. Shear and compaction band formation on an elliptic yield cap // J. Geophys. Res. - 2004. - V. 109. - B03402. - doi:10.1029/ 2003JB002633.
12. Carrol M.M. A Critical State Plasticity Theory for Porous Reservoir Rock. - Am. Soc. Mech. Eng., 1991 (117, Book G00617-1991).
13. DiMaggio F.L., Sandler I.S. Material models for granular soils // ASCE J. Engng Mech. - 1971. - V. 97. - P. 935-950.
14. de Borst R., Groen A.E. Computational Strategies for Standard Soil Plasticity Models // Modeling in Geomechanics / Ed. by M. Zaman, G. Gioda, J. Booker. - New York: John Wiley & Sons., 2000. -P. 23-50.
15. Swan C.C., Seo Y-K. Limit state analysis of earthen slopes using dual continuum/FEM approaches // Int. J. Numer. Anal. Meth. Geomech. -1999. - V. 23. - P. 1359-1371.
16. Issen K.A., Rudnicki J.W. Conditions for compaction bands in porous rock // J. Geophys. Res. - 2000. - V. 105. - No. 21. - P. 529-536.
17. Issen K.A., Challa V. Influence of the Intermediate Principal Stress on Compaction Localization Conditions // 41st US Symposium on Rock Mechanics, Colorado, June 17-21, 2006.
18. Issen K.A., Challa V. Conditions for Dilation Band Formation in Granular Materials // Proc. 16th ASCE Engineering Mechanics Conference July 16-18, 2003, University of Washington, Seattle.
19. Vajdova V, Wong T.-f., FarrellD.E., Issen K.A., Challa V Experimental Observation and Numerical Simulation of Initiation and Propagation of Compaction Bands in a Sandstone // Proc. 16th ASCE Engineering Mechanics Conference, July 16-18, 2003.
20. Katsman R., Aharonov E., Scher H. A numerical study on localized volume reduction in elastic media: Some insights on the mechanics of
anticracks // J. Geoph. Res. - 2006. - V. 111. - No. B03204. -doi:10.1029/2004JB003607.
21. Katsman R., Aharonov E. A study of compaction bands originating from cracks, notches, and compacted defects // J. Struct. Geol. -2006.- V. 28. - P. 508-518.
22. Стефанов Ю.П. Локализация деформации и разрушение в геоматериалах. Численное моделирование // Физ. мезомех. - 2002. -Т. 5 . - №5.- С. 107-118.
23. Стефанов Ю.П. Некоторые особенности численного моделирования поведения упруго-хрупкопластичных материалов // Физ. мезомех. - 2005. - Т. 8. - № 3. - С. 129-142.
24. Stefanov Yu.P. Numerical investigation of deformation localization and crack formation in elastic brittle-plastic materials // Int. J. Fract. -
2004. - V. 128(1). - P. 345-352.
25. Уилкинс М.Л. Расчет упругопластических течений // Вычислительные методы в гидродинамике. - М.: Мир, 1967. - С. 212-263.
26. Друккер Д., Прагер В. Механика грунтов и пластический анализ или предельное проектирование // Механика. Новое в зарубежной науке. Вып. 2. Определяющие законы механики грунтов. - М.: Мир, 1975. - С. 166-177.
27. Николаевский В.Н. Механические свойства грунтов и теория пластичности // Механика твердых деформируемых тел. Том 6. Итоги науки и техники. - М.: ВИНИТИ АН СССР, 1972. - С. 5-85.
28. Гарагаш И.А., Николаевский В.Н. Неассоциированные законы течения и локализации пластической деформации // Успехи механики. - 1989. - Т. 12. - № 1. - С. 131-183.
29. Николаевский В.Н. Геомеханика и флюидодинамика. - М.: Недра, 1996. - 448 с.
Поступила в редакцию 25.10.2006 г.
