CC BY
106
Некоторые особенности численного моделирования поведения упруго-хрупкопластичных материалов
Ю.П. Стефанов
Институт физики прочности и материаловедения СО РАН, Томск, 634021, Россия
В работе изложена процедура численного описания поведения упруго-хрупкопластичных материалов. Отмечены основные трудности, возникающие при численном моделировании поведения таких материалов под действием нагрузок, и методы их решения. Кратко описаны некоторые алгоритмы, как известные в литературе, так и разработанные автором, с использованием которых был решен ряд задач. Представлены некоторые результаты расчетов локализации деформации и разрушения геоматериалов.
On some features of numerical simulation of the behavior of elastic-brittle-plastic materials
Yu.P. Stefanov
Institute of Strength Physics and Materials Science SB RAS, Tomsk, 634021, Russia
The paper outlines the procedure of numerical description of the elastic-brittle-plastic material behavior. The main difficulties arising in numerical simulation of the elastic-brittle-plastic material behavior under loading and methods of their overcoming are discussed. Some algorithms known in the literature and developed by the author, which have been used to solve a number of problems, are briefly reviewed. The calculation results are given.
1. Введение
Разрушение геологических материалов связано, в первую очередь, с наличием пор и трещин различных масштабов. При любом виде нагрузки напряженно-деформированное состояние в таких материалах неоднородно. Вблизи пор и трещин могут возникнуть растягивающие напряжения, что приводит к росту трещин. Однако в силу огромного числа частиц, пор, трещин и различия масштабов возникает необходимость введения некоторых эффективных характеристик. Именно эта проблема остается одной из наиболее сложных, т.к. введение таких характеристик означает пренебрежение внутренней структурой более мелкого масштаба, которую невозможно полноценно, а зачастую и непротиворечиво учесть таким способом.
При изучении поведения хрупкопластичных материалов традиционно развиваются два направления. В первом анализируется напряженно-деформированное состояние вокруг трещин различной ориентации и рас-
сматриваются условия их роста, взаимодействия и слияния [1-4]. Во втором подходе поведение материала рассматривается как пластическое и формулируются соответствующие модели, например хрупкопластичного поведения [5-21].
В настоящий момент наиболее распространенным подходом при построении моделей роста трещин при сжатии хрупких материалов является представление о росте трещины под действием расклинивающей силы. В условиях неравноосного сжатия упругого тела, содержащего трещину, происходит смещение поверхностей трещины, в результате чего вокруг вершин возникают антисимметрично расположенные зоны растяжения и сжатия. Считается, что рост трещины такого типа носит характер отрыва и соответственно контролируется критическим коэффициентом интенсивности напряжений. Как теоретические оценки, так и анализ экспериментальных данных показывают, что трещина вытягивается в направлении оси сжатия или ортогонально оси
© Стефанов Ю.П., 2005
действия наибольших растягивающих напряжений. Одной из основных моделей роста одиночной трещины в хрупком теле, является модель Файрхерста-Кука [22]. На ее основе в ряде работ, например [23], были получены выражения для длины такой трещины, ее раскрытия, коэффициента интенсивности напряжений. Однако аналитические выражения позволяют провести оценку состояния трещины в хрупком материале для строго оговоренных условий. В случае близости границы раздела разнородных материалов эти условия становятся достаточно сложными. Кроме того, одной из основных проблем, связанных с исследованием роста трещин, является то, что все такие оценки справедливы лишь до начала роста трещины. С момента страгивания трещины (даже без учета динамических явлений), строго говоря, мы имеем дело с новым объектом. Поскольку геометрия тела изменилась, необходимо заново решать задачу о напряженном состоянии уже с новыми граничными условиями, что практически невозможно осуществить аналитическими методами. В связи с этим наиболее подходящий путь решения таких задач — это применение численных методов, которые позволяют рассчитать напряженно-деформированное состояние в ходе развития трещины.
В том случае, если число трещин велико, а напряженное состояние препятствует их раскрытию и формированию магистральных трещин, т.е. до тех пор, пока в целом образец сохраняет сплошность, наиболее подходящий способ описания поведения — это использование соответствующей модели упруго-хрупкопластического поведения. Независимо от условий нагружения и проявления видимого макроскопического поведения прочность геологических материалов традиционно определяется двумя параметрами: когезией, или сцеплением, и трением частиц между собой и на поверхностях трещин. Чтобы наиболее адекватно описать поведение материала при различных условиях нагружения, модель должна учитывать накопление повреждений и связанное с ними макроскопическое изменение объема.
Наибольшее распространение получили различные варианты обобщения модели Друккера-Прагера с неассоциированным законом течения. Широко используются модели Кулона-Мора, Хьюка-Брауна, а также модели, использующие различные формы третьего инварианта напряжений, например Лоде, Лоде-Дункана. Описание ряда таких моделей можно найти в работах [20, 21].
Для описания процесса деформирования будем использовать комплексный подход, согласно которому считаем, что в ходе неупругой деформации происходит накопление повреждений, снижающих прочность материала, а раскрытие трещин происходит под действием растягивающих напряжений. Описание неупругой деформации, при сохранении сплошности, будем осуществлять, используя подход и математический аппарат
теории пластического течения. Поэтому в дальнейшем под пластической деформацией будем понимать неупругое поведение, независимо от его природы. Учет образования и роста трещин при численном моделировании будем осуществлять явным образом, с формированием свободных поверхностей при помощи специального алгоритма разделения узлов сетки. Такой способ описания трещинообразования обеспечивает автоматический учет концентрации напряжений в вершинах трещин, независимо от их числа.
2. Определяющие соотношения для описания поведения хрупкопластичных материалов
Будем считать, что полные деформации е- состоят из упругих е® и пластических ер-. Аналогичное разложение примем и для скорости деформации:
еУ- = ее+Ц. а)
Связь между напряжениями и деформациями, до наступления пластической деформации, опишем гипоуп-ругим законом:
= - Р8- + ,
—- - 2ц[ єгу — гекк8„ Ш I 1 3 1
Ш-
(2)
(3)
(4)
(5)
Р = - К^.
V
Здесь а- — компоненты тензора напряжений; я- — компоненты девиатора тензора напряжений; Р — среднее давление; г- — компоненты тензора деформаций Коши:
(6)
иг — компоненты вектора перемещения; ю- — компоненты тензора скоростей вращения:
= \(иі,- -
(7)
К и ц — модули сжатия и сдвига соответственно.
Неупругая деформация определяется в соответствии с заданными поверхностью нагружения и законом течения:
/е, х) = 0, (8)
Эя
(9)
где f— поверхность текучести; х — параметр упрочнения; g — пластический потенциал; X — множитель, определяемый в ходе процесса деформации.
Чем сложнее функция текучести (8), тем труднее определить значения параметров и тем сложнее должны быть эксперименты для их определения. Поэтому целесообразно использовать наиболее простой вид функции
текучести, который бы позволил с приемлемой точностью описать свойства материала в интересующих условиях.
Наиболее простой формой уравнения, задающего поверхность текучести для чувствительных к давлению хрупких материалов, являются зависимости, включающие первый инвариант напряжений J1 [24]:
/ = а11 + Jl/2 - У, (10)
f = аYJl + J2 - Y2.
(11)
Здесь У — предел текучести при сдвиге; а — некоторый
коэффициент; J1 = а1 + а2 + а3 = -3Р; J 2 = — (.—■,■).
2
Первое из этих уравнений является обобщением условия Кулона-Мора и определяет конус Мизеса-Шлей-хера. Второе уравнение задает параболоид вращения. Условие такого типа использовано, например, в работах
[11, 25].
Наиболее распространенным видом уравнения для описания поведения геоматериалов за пределом упругости остаются различные формы записи уравнения Кулона-Мора |т| = С + tgф ап, где ф — угол внутреннего трения; т — наибольшее касательное и ап — нормальное напряжения; С — коэффициент сцепления, а также его обобщения с неассоциированным законом течения [26, 27]:
f = аі - аз^ + 2СдРф ,
где Nф = (1 + sin ф)/(1 - sin ф),
g = а1 -аз Nv,
(12) (13)
где = (1 + sin у)/(1 - sin у); у — угол дилатансии.
В работе Друккера и Прагера [5] рассмотрено применение ассоциированного закона течения для обобщенного условия Кулона-Мора в форме Мизеса-Шлейхера: aJ + J 22 = Y, (14)
где a — коэффициент внутреннего трения; Y — сдвиговая прочность материала, или, по аналогии с условием Кулона-Мора, сцепление.
При использовании данной модели, ставшей одной из наиболее распространенных, приращения пластической деформации определяются выражением:
•Л + Jt
(15)
Таким образом, пластическая деформация имеет не только сдвиговой, но и объемный характер, между их характеристиками возникает связь: = ба/р12. Усло-
вие такого типа было рассмотрено Новожиловым в [28].
Однако при использовании ассоциированного закона течения не удается получить количественное согласование с наблюдаемыми данными при деформации геоматериалов. Величина объемной пластической деформации оказывается существенно завышенной, т.к. коэффициент дилатансии Л, связывающий объемную и сдви-
говую части деформации в виде соотношения 11р = = 2Л1р12, равен 3 а.
Введение второго независимого параметра — коэффициента дилатансии — осуществляется при обобщении модели [13, 19, 20], с введением пластического потенциала в форме:
я = и+J2/2. (16)
В этом случае приращения пластической деформации определяются выражением:
, * )
dep- = dA
Ч
р8
У + 2jf
(17)
В данной работе использована модель Николаевского с независимым соотношением между объемной и сдвиговой составляющими пластической деформации [8, 9]. В качестве условия пластичности применяется условие Мизеса-Шлейхера в форме:
J + J f = Y.
(18)
Пластический потенциал записывается в виде:
Л
а
g(а-) = J2 + — J 2Y-— J1 l + const.
.у/ „,.з„ ^ (19)
Приращения пластической деформации определяются выражением:
2
S і +— Л
гі З
y-а j
з
8г,
dA.
(20)
Связь между объемной и сдвиговой пластической деформациями определяется в виде:
1р = 2Л72р1/2 (21)
где Л — коэффициент или скорость дилатансии.
3. Разупрочнение и разрушение
Для численного моделирования разделим процессы локализации деформации и разрушения без видимого нарушения сплошности и формирование явных трещин. Эти процессы будем описывать двумя различными подходами, используя параметры, отражающие накопление повреждений.
Будем считать, что сдвиговая прочность материала зависит от накопленной пластической деформации:
У = Уо (1 + А(е) - Ве (е)), (22)
где А(е) — функция, описывающая упрочнение материала; Бе (е) описывает накопление повреждений и разупрочнение, е = 21р1/2.
Примем линейную зависимость для слагаемого, характеризующего упрочнение,
А(е) = И Ц (23)
е
и квадратичную — для разупрочнения по мере накопления повреждений, проявляющихся в виде дилатансии материала,
(
Ое (е) = 2И
/р/Л
= 2И
е
2
=А
е' ( е *
е
(24)
где к — коэффициент упрочнения; е* — критическая деформация, после достижения которой преобладает деградация материала.
При действии растягивающих напряжений ап > 0 возможно образование и рост трещин. Условие разрушения запишем в виде:
аей- = ап (1 - а) + ата < а*(1 - Бе) для ап > 0, (25) где ап, ат — нормальная и тангенциальная компоненты вектора напряжений; 0 < а < 1 — параметр, описывающий вклад нормальных и касательных напряжений;
% ч
а — предельное значение напряжений.
Таким образом, принимаем, что разрушение и раскрытие трещин происходят при наличии и под действием растягивающих напряжений, которые присутствуют в локальных областях неоднородной среды практически при любом виде нагрузок. В ходе пластической деформации происходит накопление повреждений, связанное с образование и ростом микротрещин и пор. Развитие и накопление этих повреждений проявляется в объемной пластической деформации и разупрочнении среды.
4. Процедура расчета
Описание процесса деформации осуществляется путем численного решения системы уравнений механики сплошной среды:
- =Ри г, (26)
Р + Ри, ■ = 0. (27)
Система уравнений решается при заданных начальных и граничных условиях. Замыкают систему уравнений модельные соотношения (1)-(7), (18)-(25) упругохрупкопластического тела, представленные в предыдущем разделе.
Задачи решались в двумерной постановке для условий плоской деформации. Для решения системы уравнений использовалась явная конечно-разностная схема из [29].
На каждом последовательном шаге по времени приращения пластической деформации будут пропорциональны разнице между напряжениями, рассчитанными по упругому закону (обозначим их *), и напряжениями, соответствующими предельной поверхности (текучести). При переходе на следующий шаг по времени выполним предварительный расчет напряжений, соответствующих приращению полных деформаций по упругому закону:
................ (28)
(29)
(*П+1)* = + (д^+1)*,
После этого проверим условие (18), т.е. проверяем, не оказалось ли напряженное состояние за предельной поверхностью, в запрещенной области.
Эту проверку удобно осуществить, записав условие в виде:
С = а Г* + J 2*12 -3 1 2
У.
(30)
„П+1\*
Здесь
Jl* = J 2 = J2(sГ Г. (31)
Если С < 0, то материал в данной ячейке расчетной сетки в текущий интервал времени находится внутри предельной поверхности, т.е. в упругом состоянии. Тогда рассчитанное напряженное состояние соответствует истинному и можно переходить к дальнейшему расчету на следующем слое по времени.
Если С > 0, то материал перешел в пластическое состояние и необходимо скорректировать напряжения в соответствии с законом течения так, чтобы они оставались на предельной поверхности, а также рассчитать приращения пластической деформации. Соотношения, необходимые для осуществления такой процедуры, были получены в [30]. Запишем их в виде:
3КЛС
Т = Т --
Jl = Jl
(32)
аКЛ + ц
J У2 = J2*1/2-----^. (33)
2 2 аКЛ + ц
Имея значения второго инварианта девиатора напряжений, вычислим компоненты напряжений как:
Т V2
* ° 2
* - = ^лг.
•Т 2
Параметр X определяется выражением:
х=^12---------------------С-.
J *12(аК 2Л + 2ц) - 2цС
(34)
(35)
(Ри+1)* = Рп + (ДР”+1)*.
Теперь по формуле (20) могут быть вычислены приращения пластической деформации.
Следует заметить, что при Л = а мы имеем модель Друккера-Прагера и при Л = а = 0 — модель Прандт-ля-Рейса со стандартной процедурой приведения напряжений на круг текучести [29].
Значения У, а, Л не остаются постоянными в ходе деформирования. Понятно, что они изменяются по мере относительного смещения поверхностей микротрещин. Используемый метод допускает выполнение расчетов при изменении их значений на каждом временном интервале вычислений, по мере изменения состояния материала. Параметры, определяющие поверхность течения и направление вектора пластической деформации, У, а и Л приобретают смысл функций от накопленной сдвиговой и/или объемной пластической деформации, значения которых постоянны в течение одного временного шага интегрирования.
2
*
е
Ац V
—,
)Т —"м г ^ ■
г* - >0
Рис. 1. Вид предельной поверхности с ограничением в области растяжения и ее изменение в ходе деформации. Угол внутреннего трения и предел текучести зависят от накопленной пластической деформации
Несмотря на то, что выражения (32)-(35) получены для идеального случая с постоянными коэффициентами У, а и Л, в процедуре расчета используются приращения напряжений и деформаций за малые промежутки времени. Поэтому ничто не мешает рассматривать данные коэффициенты как функции накопленной пластической деформации, считая их постоянными в течение шага по времени.
Поверхность, задаваемая уравнением (18), соответствует конусу Кулона-Мора. Если принять, что ее параметры есть функции накопленной пластической деформации (22), то в ходе деформации поверхность будет меняться. Вид этой поверхности в координатах (Р, 312) и ее изменение в ходе деформации показаны на рис. 1. Принимая в качестве параметра состояния накопленную пластическую (неупругую) деформацию, можно задать величины угла внутреннего трения ф и предела текучести (сдвиговой прочности) У как функцию этого параметра. Тогда поверхность текучести будет меняться в ходе развития деформации. Она приобретает смысл мгновенной поверхности, соответствующей данному уровню накопленной неупругой деформации, и может быть различной в различных точках среды.
5. Описание роста трещин
Моделирование процессов деформации с учетом разрушения материала относится к разряду наиболее
алгоритмически сложных задач. Для описания образо-
вания и роста несплошностей необходимо применение
специальных процедур, которые позволяли бы не только
определить время и место возникновения или роста тре-
щины, но также учитывать образование и наличие но-
вых поверхностей внутри расчетной области в ходе дальнейшего расчета. Для определения местоположения и момента образования несплошностей применяются различные модели и критерии разрушения. Обзор критериев разрушения для случая динамического нагружения тел можно найти в работе [31].
Для расчета геометрии и напряженно-деформированного состояния в разрушенной области используют два различных подхода. Первый из них — неявный, состоит в коррекции напряжений и/или даже удалении расчетных ячеек в зоне разрушения. Здесь необходимо отметить одну из основополагающих работ [32], в которой принимаются во внимание ориентация и раскрытие образующихся трещин. Второй — явный, предполагает введение свободных поверхностей, соответствующих берегам трещины. В этом случае изменения напряженно-деформированного состояния, связанные с разрушением, учитываются автоматически. Существуют различные схемы описания роста трещины, как путем перестройки сетки вблизи вершины, например [33], так и по границам расчетных ячеек [34-36]. Причем используются различные процедуры расчета — от скачкообразного до плавного перемещения. Обзор численных подходов описания распространения трещин можно найти в [37-39].
Будем считать, что процессы разрушения локализованы по границам ячеек. Чтобы описать образование и рост трещин, используется процедура разделения узлов расчетной сетки [34-36].
Новые поверхности трещин образуются по границам расчетных ячеек (рис. 2). Проверка условия разрушения осуществляется во всей расчетной области на границах между ячейками. При выполнении критерия происходит раскрытие «элементарной» трещины, в качестве которой выступают границы соответствующих расчетных ячеек. Поскольку при решении используется прямоугольная сетка, каждая возможная элементарная трещина может быть ориентирована только в двух взаимно перпендикулярных (пока сетка не деформирована) направлениях. На рис. 2 показаны все возможные конфигурации элементарных трещин. Таким образом, автоматически определяются не только местоположение разрушаемой области, но и ориентация элементар-
Рис. 2. Возможные конфигурации трещин при разделении узлов расчетной сетки
Рис. 3. Нежелательное искажение формы ячеек расчетной сетки в виде «песочных часов» (а), дополнительное разбиение ячеек на треугольники в «треугольной» вязкости (б), углы ю, изменение которых рассматривается при использовании «угловой» (в) и «моментной» (г) вязкостей
ных трещин, рост и взаимодействие которых формируют трещины большего масштаба.
Для выполнения такой операции каждый узел сетки заменяется группой из четырех узлов. До тех пор пока материал остается сплошным, узлы в группе объединены, а схема расчета сохраняется прежней. В случае необходимости группа распадается, и мы имеем не один узел, а два, три или четыре узла, в соответствии с конфигурацией трещин в разрушенной области. При этом не меняется схема расчета параметров напряженно-деформированного состояния, т.к. они вычисляются в ячейках сетки.
После разделения узел расчетной сетки оказывается принадлежащим образованной поверхности. Для дальнейшего расчета движения необходимо задание на всех созданных таким образом границах соответствующих условий. В простейшем случае, при раскрытии трещины, это условия свободной поверхности. В случае сближения берегов, чтобы избежать перехлеста лежащих по разные стороны от трещины ячеек, необходимы проверка и применение контактных условий. Простейшим вариантом таких условий может быть совмещение узлов. Эти условия могут быть использованы в случае трещины отрыва, если относительное поперечное смещение невелико. В противном случае возникает необходимость применения более сложных контактных условий. Таким образом, после разделения узлов на вновь образованных поверхностях задаются условия свободной поверхности и условия взаимного непроникновения.
Расчет взаимодействия контактных границ чаще всего осуществляется по схеме коррекции движения узлов, предотвращающей взаимное проникновение. Проблема расчета взаимодействия сеток оказывается одной из наиболее сложных при исследовании контактных задач и разрушения. Реализация алгоритмов взаимодействия существенно усложняет численный код и значительно увеличивает временные затраты на вычисления. Обзор
наиболее распространенных алгоритмов и подходов к решению данной задачи можно найти в работе [40].
Реализация граничных условий на контактирующих поверхностях трещины осуществлялась по схеме коррекции движения узлов расчетных ячеек. Используемый алгоритм реализации контактных условий аналогичен алгоритму, описанному в [33, 41]. Некоторые особенности и возможности использованных алгоритмов разделения узлов и контактного взаимодействия поверхностей трещины описаны в [30, 42-46]. Хорошие результаты были получены также с использованием алгоритма, изложенного в [47, 48], в котором расчет взаимодействия осуществляется на этапе расчета скоростей, учитывая силы, действующие по разные стороны от границы. Данный алгоритм был использован для реализации граничных условий на границе раздела расчетных областей, описываемых конечно-разностным и дискретным методами.
Еще одной проблемой, возникающей при моделировании, является проявление численных эффектов, не соответствующих физике процесса. Например, хорошо известно, что треугольные ячейки имеют избыточную жесткость, что может исказить результаты при больших деформациях. В противоположность этому «жесткость» четырехугольных ячеек недостаточна. Одной из самых серьезных проблем применения методов с прямоугольными ячейками являются так называемые «песочные часы» — трапецевидное изменение формы ячеек (рис. 3, а), которое не приводит изменению напряженного состояния, т.к. вычисленное значение деформации не меняется. Такой тип деформации может очень быстро привести к невозможности продолжения расчета. Для предотвращения такого типа деформации используются «угловые» [32] и «треугольные» [49] формы искусственной вязкости.
Применение треугольного типа вязкости фактически означает дополнительное разбиение ячеек на треуголь-
ники, окружающие каждый узел (рис. 3, б). В этих ячейках рассчитываются скорость деформации и значения искусственной вязкости — дополнительных «вязких» напряжений, действующих на узлы сетки [49]. Искусственная вязкость такого типа хорошо работает при расчете больших деформаций, однако приводит к достаточно ощутимым затратам на вычисления.
При использовании искусственных вязкостей углового типа [32] (рис. 3, в) принимается во внимание разница изменения углов противоположных граней. Ее вычисление также требует существенных вычислительных затрат. Значительно экономичнее оказывается алгоритм «моментной» вязкости, в котором рассматривается изменение отношения проекции длины противоположных граней ячейки (рис. 3, г). Такой алгоритм особенно эффективен при расчете малых деформаций, например, задач об излучении и распространении упругих волн.
Однако не существует универсальных рецептов по выбору того или иного типа вязкости и значений их коэффициентов. Для каждого типа задач приходится выбирать наиболее удобный вариант. Например, для задач, в которых имеют место большие значения деформаций и повороты ячеек, хорошо работает алгоритм треугольной вязкости. Для задач, в которых деформации ячеек невелики, вероятно, наиболее подходящими будут угловые типы вязкости. Здесь также хорошо работают алгоритмы, где рассматривается не поворот, а отношение длин противоположных граней.
В числе проблем, которые возникают при использовании схемы [29], следует выделить осцилляции, связанные с точечным распределением масс и порядком аппроксимации уравнений. Для устранения данных эффектов хорошо работают известные типы искусственных вязкостей: объемные квадратичные и линейные, а также тензорные или вязкости типа Навье-Стокса [29]. Наилучшие результаты получаются при использовании комбинации вязкостей различного типа. В работе [50] для некоторых случаев предлагается также использовать алгоритмы сглаживания (размазывания), а также фильтрации высокочастотных колебаний при анализе результатов.
6. Обсуждение результатов
На примерах решения ряда модельных задач рассмотрим некоторые особенности процесса деформирования при нагружении тел, содержащих несплошности.
Пусть тело, содержащее трещину, подвергается воздействию ударной нагрузки. Движение поверхностей трещины при нормальном падении волны растяжения, а также сопутствующие явления дифракции и излучения упругих волн вершиной подробно рассмотрены в [46]. Поскольку в этом случае берега трещины расходятся и не контактируют друг с другом, очевидно, что для расчета не требуется учитывать их взаимодействие. То же
А ” 1 В
Рис. 4. Схема нагружения образца с трещиной
самое можно сказать и в случае произвольного угла падения волн растяжения на поверхности трещин. В связи с этим имеет смысл рассмотреть задачи о падении на трещины волн сжатия и сдвига, когда взаимодействие поверхностей играет существенное значение [42].
6.1. Трещина под действием комбинированной нагрузки
Рассмотрим задачу о поведении тела под действием комплексной ударной нагрузки сжатия и сдвига. На рис. 4 показаны геометрия образца и схема его нагружения. Граничные условия имеют вид: их (х, у, 0 = /х (/), иу (х, у, г) = /у (4 (х, у) е АВ.
На границах AD и ВС задавались условия симметрии, на DC — условие свободной границы.
При таком воздействии мы будем наблюдать распространение двух типов волн: продольной (волны сжатия) и поперечной (волны сдвига), фронты которых на некотором расстоянии от границы АВ разделятся. Оба фронта являются плоскими и распространяются параллельно плоскости трещины. После достижения поверхности трещины, волна сжатия пройдет сквозь нее без каких-либо искажений, т.е. в этом случае материал в окрестности такого разреза должен вести себя как сплошной. Это значит, что до прихода сдвиговой волны точки, принадлежащие противоположным берегам разреза, двигаются совместно. Наличие трещины обнаружит только волна сдвига. Для этой волны, с фронтом параллельным плоскости разреза, при отсутствии трения, поверхности трещины являются свободными. Обратившись к рис. 5, мы видим именно такое движение. Различие в скоростях движения возникает с приходом поперечной волны. С этого момента у одного из берегов появляется отличная от нуля х-компонента скорости. Несколько позже возникает и ^-компонента, но уже на обеих поверхностях трещины. Появление смещений в этом направлении связано с приходом возникших в вершинах дифрагированных волн.
Рис. 5. Компоненты скорости точек, лежащих на противоположных берегах трещины, при последовательном падении волн сжатия и сдвига
В области трещины волна сдвига отражается, не меняя своего типа и направления смещения частиц. В то же время, вершины трещины становятся источником искривления фронта и излучения волн различных типов. Свободное смещение (сдвиг) поверхностей трещины приводит к возникновению антисимметрично расположенных зон растяжения и сжатия. В вершинах происходит концентрация напряжений. На рис. 6 приведены картины распределения компонент напряжений. На них хорошо видно, что в соответствии со сказанным выше, в области, которой волна сдвига еще не достигла, напряженное состояние является однородным. Возмущения, связанные с взаимодействием продольной волны сжа-
тия и трещины, отсутствуют, что соответствует физике процесса.
6.2. Наклонное падение волны сжатия на поверхность трещины
При наклонном падении волны сжатия взаимодействие и скольжение поверхностей происходят по мере достижения фронтом волны поверхности трещины. Точка пересечения фронта с трещиной становится движущимся источником излучения упругих волн. Причем очевидно, что эта кажущаяся скорость — скорость источника — будет превышать скорость самой волны: V = = Vplcos 0, где 0 — угол наклона трещины относительно направления распространения фронта волны нагрузки. Движение источника приводит к образованию конических волн продольного и поперечного типа. На рис. 7 показаны возмущения поля скоростей при наклонном падении плоской волны сжатия на трещину в последовательные моменты времени. При данной геометрии фронт головной волны продольного типа образует угол 45° к поверхности трещины. Поэтому позади трещины он лежит на линии фронта волны нагрузки и сливается с ним. Фронт поперечной конической волны, ввиду меньшей скорости ее распространения, образует более острый угол (=36°) к поверхности трещины в соответствии с законом Снеллиуса. Для иллюстрации описанного выше процесса на рис. 8 приведено распределение х- и у-компонент напряжений. Также, как и в рассмотренных выше случаях, к концентрации напряжений приводит
сжатия и сдвига в образце с трещиной в области, отмеченной пункти-
■: ЛЛ1 - ■ Г---1Г ■. л/-" -т.- ром на рис. 4
Рис. 7. Возмущение поля скоростей при взаимодействии плоской волны сжатия с наклонной трещиной: W — фронт падающей волны, Р -продольные и S — поперечные волны
сдвиг поверхностей трещины. Рост напряжений в вершине начинается после того, как ее достигнет волна нагрузки.
6.3. Деформирование тела с множеством скользящих слоев
Для иллюстрации работы расчетных алгоритмов приведем результаты моделирования поведения тела, содержащего множество свободных поверхностей. В этом случае взаимодействие поверхностей полностью определяет картину деформирования.
Рассмотрим поведение тела, содержащего множество вложенных Г-образных слоев, в условиях медленного (квазистатического) сжатия. Поведение материала описывалось в рамках упруго-идеальнопластической модели. Геометрия образца и схема его нагружения показаны на рис. 9, а. На рис. 9, б, где приведена полученная деформированная сетка, хорошо видны смещения слоев друг относительно друга и возникающие при этом пустоты. Часть слоев изгибается, при этом удается избежать проблем, связанных с деформированием ячеек типа «песочные часы». Легко можно увидеть, что возникающие
в процессе деформирования внутренние пустоты приводят к увеличению общего объема тела.
6.4. Разрушение упруго-хрупкопластичного образца с включением
Наличие включений приводит к сложному напряженному состоянию. Вблизи пор и трещин возникают растягивающие напряжения, что приводит к росту трещин отрыва. Данный процесс можно проиллюстрировать на примере однородного образца, содержащего включение более «мягкого» материала или отверстие в центральной части. В этом случае локализация деформации и разрушение инициируются в этой области. Например, на рис. 10 представлены численные расчеты разрушения таких образцов с явным учетом образования трещин. Условия нагружения соответствовали вертикальному сжатию вдоль оси образца с идеальным скольжением на торцах. Развитие трещин инициировалось в области отверстия в центре образца.
При деформировании среды, содержащей включения более жесткого материала, неоднородность напряженно-деформированного состояния приводит к пово-
Рис. 8. Распределение компонент напряжений при взаимодействии плоской волны сжатия с наклонной трещиной
'/7777777777777777777777/,
Рис. 9. Схема нагружения образца, содержащего множество Г-образных слоев (а) и рассчитанная картина деформирования (б)
роту включения, который сопровождается возникновением зон растяжения, и образованию трещин. На рис. 11 показан фрагмент расчетной сетки в области разрушения для такого случая.
6.5. Вдавливание жесткого штампа
На рис. 12 представлены результаты численного моделирования процесса деформации упругопластического слоя, лежащего на жестком основании, при вдавливании жесткого штампа. В зависимости от свойств материала, аналогичные условия нагружения приводят к принципиально различным картинам деформирования.
Рис. 10. Картина разрушения образцов с центральным отверстием при сжатии
Рис. 11. Поворот жесткого включения и разрушение вмещающей среды в результате сдвига
Рис. 12. Распределение интенсивности деформации вблизи области приложения нагрузки при вдавливании жесткого штампа в слой среды, лежащий на жестком основании, при различных значениях внутреннего трения и дилатансии: а = Л = 0 (а); а = 0.25, Л = 0.15 (б); а = 0.546, Л = 0.12 (в)
Например, рис. 12, а получен с использованием модели Прандтля-Рейса (с условием текучести Мизеса). Это наиболее распространенная модель, дающая хорошие результаты при описании пластичных материалов. Однако она не учитывает влияние давления на прочность среды и объемные изменения, связанные с пластической деформацией. Для описания неупругого поведения геологических сред более адекватные результаты позволяют получить модели, учитывающие внутреннее трение и дилатансию. При малых значениях этих параметров картина деформирования остается похожей на первую (рис. 12, б). Вариация коэффициентов (что соответствует рассмотрению различных материалов) показывает, что при определенных значениях коэффициентов среда становится склонной к локализованному развитию деформации. Возможно возникновение системы полос локализованной деформации, распространяющихся вглубь слоя (рис. 12, в), что соответствует более хрупкому поведению среды.
Таким образом, в зависимости от выбора модели и значений параметров, т.е. учета тех или иных особенностей поведения материала, могут быть получены различные картины деформирования. Это остается достаточно серьезной проблемой при описании поведения геологической среды, тем более что свойства геологических материалов остаются мало изученными и их поведение сильно зависит от условий нагружения. Отмеченные особенности решений, полученных с использованием различных моделей, согласуются с результатами и выводами работы [51].
6.6. Локализация деформации и разрушение
при сжатии образцов упруго-хрупкопластичного
материала
Рассмотрим поведение однородного образца, механические характеристики которого соответствуют макро-
Об
■04
1
' % ■ С-.2
оа
О Q.2 04 й.6 ОВ
if*, Vi
Рис. 13. Параметры модели как функции накопленной пластической деформации. Кривые соответствуют состоянию ячейки в центре образца
скопическим свойствам песчаника [52]: р0 = 2.2 г/см3, K = 12.28 ГПа, ц = 5.346ГПа, 70 = 9.04МПа, а0 = 0.546, Л 0 = 0.48 в условиях вертикального сжатия вдоль оси образца. Зависимости параметров Y, а и Л от неупругой деформации показаны на рис. 13.
При сжатии прямоугольного однородного образца на величину осевой деформации до разрушения и на картину разрушения существенно влияют условия на торцевых поверхностях. С увеличением коэффициента трения на торцах локализация деформации и разрушение происходят при меньшей степени общей осевой деформации и соответственно при меньшей величине нагрузки (рис. 14). Несколько отличается и картина полос локализации деформации.
Развитие деформации в условиях идеального скольжения на торцах до достижения предельного уровня протекает однородно. При достижении предельного уровня деформации начинают формироваться полосы локализации (рис. 15, а, б).
0 09
000
ш
^ 004
t
ОН
o.i о г oj од AV7V, к
Рис. 14. Кривые зависимости напряжений а = а! — ©2 от осевой (а) и объемной (б) деформации образца песчаника при различных условиях на торцевых поверхностях. Кривые 1 соответствуют расчетам для образцов с «фиксированными» (трение велико) торцами, кривые 2, 3 — при идеальном скольжении. Кривая 4 соответствуют экспериментальным данным работы [52]
Рис. 15. Формирование сетки полос (а, б) и типичная картина локализации деформации (в) в образце песчаника при обжимающем давлении ©2 = 5 МПа. Однородный образец с идеальным скольжением на торцах
На начальном этапе в образце образуется сетка полос локализованной деформации. На этом этапе места зарождения полос локализации деформации определяются малыми по величине осцилляциями напряжений и деформаций, связанными с волновым фоном и динамической формулировкой задачи, а также ошибкой вычислений. В выполненных расчетах их величина не превышает 0.01 % от максимального уровня действующих в образце напряжений. По мере дальнейшего деформирования развитие локализации происходит строго определенным образом. Развиваются лишь полосы, обусловленные наличием особенностей структуры и геометрией образца. Поскольку в рассмотренном случае моделировалось поведение однородного образца, полосы локализованной деформации развивались лишь от угловых точек. Например, на рис. 15 показано распределение интенсивности деформации в однородном образце на различных этапах развития локализации при граничных условиях, соответствующих идеальному скольжению на торцевых поверхностях.
Увеличение трения на торцевых поверхностях приводит к более интенсивному развитию деформации в углах образца и, как следствие, к более быстрому появлению полос локализованной деформации. В этом случае локализация деформации и разупрочнение происходят на более ранней степени деформирования. Явный учет образования трещин не приводит к существенному изменению результатов, т.к. трещины формируются в зоне локализованной деформации.
Одной из особенностей поведения такого рода материалов является эффект дилатансионного упрочнения. В работах [7, 9] показано, что локализация дефор-
мации возможна при определенных значениях положительного коэффициента упрочнения. В то же время, расчеты показывают, что эффект упрочнения может наблюдаться при нулевом и даже отрицательном коэффициенте деформационного упрочнения. Например, на рис. 13 видно, что а = а1 — а2 продолжает свой рост до достижения осевой деформации е = 0.6 %, тогда как максимальное значение У было достигнуто при е = 0.165% и далее снижалось. Данный факт объясняется высокой скоростью дилатансии и соответствующим ростом давления.
Использование рассмотренных подходов для описания поведения образцов горной породы как идеальной однородной среды позволяет получить хорошее совпадение картин разрушения, а также отдельных зависимостей напряжений и деформаций с экспериментальными наблюдениями. Однако необходимо заметить, что некоторые особенности поведения таких материалов в рамках однородной среды описать не удалось. Более того, результаты расчетов противоречат некоторым экспериментальным наблюдениям.
На поведение горных пород большое влияние оказывают величина бокового обжатия и условия на торцевых поверхностях [10]. Боковое обжатие образца ограничивает возможности роста трещин, так что при определенных условиях материал может вести себя как пластичное тело. С увеличением трения на торцевых поверхностях породы показывают большую прочность, и наоборот, наличие смазки, обеспечивающей скольжение на торцах, приводит к снижению прочности и нередко к раскалыванию образца. Имеются два традиционных объяснения этих особенностей: трение приводит к появ-
Рис. 16.Кривые зависимости напряжений а = ©! —©2 от осевой деформации неоднородного образца при различные условиях на торце-выгс поверхностях. Кривые 1, 2 соответствуют расчетам для образцов с «фиксированными» (трение велико) торцами и при условии идеального скольжения соответственно
лению сжимающих напряжении, которые препятствуют прорастанию трещин к торцам, аналогично обжимаю-щеИ нагрузке; смазка поступает в поры образца, что вызывает ряд дополнительных эффектов, на чем в данноИ работе мы останавливаться не будем.
При рассмотрении поведения однородных образцов в условиях сжатия растягивающие напряжения не возникают (при малоИ скорости нагружения, когда динамические эффекты не проявляются). В этом случае очагами локализации деформации остаются только особенности геометрии, а также условия приложения нагрузки.
В случае неоднородного пористого образца начинает эффективно проявляться отрывноИ механизм разрушения. При отсутствии трения на торцах образца развитие неупругоИ деформации начиналось несколько раньше. Однако в дальнеИшем в этом случае пластическая деформация имела более рассеянныИ характер, чем в образце с фиксированными торцами (большим трением), где пластическая деформация быстрее развивалась в центральноИ части. Наибольшее увеличение «прочности» при наличии трения было получено в образцах, имеющих большиИ размер пор и трещин вблизи торцевых поверхностеИ. Прочность образцов с фиксированными торцами оказывалась выше на 10-15 % (рис. 16, 17).
7. Заключение
В работе представлен краткиИ обзор проблем, связанных с численным моделированием поведения упруго-хрупкопластичных материалов под деИствием нагрузок, и методов их решения. Приведены некоторые алгоритмы, как известные в литературе, так и разработанные автором, с использованием которых был решен ряд задач.
Использование моделеИ упруго-хрупкопластических материалов позволяет описать многие особенности их
Рис. 17. Картины локализации деформации в неоднородном пористом образце при идеальном скольжении на торцевых поверхностях (а) и «фиксации» торцов (б)
поведения. Такие расчеты могут быть полезны для проверки адекватности модели, а также для объяснения некоторых особенностей поведения. В то же время, использование упруго-хрупкопластичной модели для описания сложных неоднородных пористых и трещиноватых сред как однородной среды не всегда позволяет учесть некоторые особенности их поведения. Для наиболее полного описания особенностей поведения неоднородных пористых и трещиноватых материалов, необходимо принимать во внимание особенности их структуры, а также рост трещин и взаимодействие их поверхностей.
Работа выполнена при поддержке РФФИ (гранты № 05-05-64659-а и № 04-05-64547-а), а также проекта 6.5.2 Программы специализированных отделений РАН.
Литература
1. Райс Дж. Механика очага землетрясения // Механика. Новое в зарубежной науке. Вып. 28. Механика очага землетрясения. - М.: Мир, 1982. - С. 10-132.
2. Gol’dstein R.V., KaptsovA.V Formation of fracture structures of weakly interacting cracks // Mechanics of Solids. - 1982. - 17(4). - P. 157166.
3. Cherepanov G.P. Mechanics of brittle fracture. - New York: McGraw Hill, 1979. - 939 p.
4. Germanovich L.N., Cherepanov G.P On some general properties of strength criteria // Int. J. Fract. - 1995. - V. 71. - P. 37-56.
5. Друккер Д., Прагер В. Механика грунтов и пластический анализ или предельное проектирование // Механика. Новое в зарубежной науке. Вып. 2. Определяющие законы механики грунтов. - М.: Мир, 1975. - С. 166-177.
6. Райс Дж.Р. Локализация пластической деформации // Теоретическая и прикладная механика: Труды XIV Межд. конгресса IUTAM, 30.08-4.09.1976, Delft, Нидерланды / Под ред. В.Т. Койте-ра. - М.: Мир, 1979. - С. 439-471.
7. Rudnicki J.W., Rice J.R. Condition for localization of plastic deformation in pressure sensitive dilatant materials // J. Mech. and Phys. Solids. - 1975. - V. 23. - No. 6. - P. 371-390.
8. Николаевский В.Н. Механические свойства грунтов и теория плас-
тичности // Механика твердых деформируемых тел. Т. 6. Итоги науки и техники. - М.: ВИНИТИ АН СССР, 1972. - С. 5-85.
9. Гарагаш И.А., Николаевский В.Н. Неассоциированные законы течения и локализации пластической деформации // Успехи механики. - 1989. - Т. 12. - № 1. - С. 131-183.
10. Поль Б. Макроскопические критерии пластического течения и хрупкого разрушения // Разрушение, Т. 2. Математические основы теории разрушения / Под ред. Г. Либовица. - М.: Мир, 1975. -С.336-520.
11. Драгон А., Мруз 3. Континуальная модель пластически хрупкого поведения скальных пород и бетона // Механика деформируемых твердых тел. Направления развития. - М.: Мир, 1983. - С. 163188.
12. Бажант 3. Эндохронная теория неупругости и инкрементальная теория пластичности // Механика деформируемых твердых тел. Направления развития. - М.: Мир, 1983. - С. 189-229.
13. Jirasek M., Bazant Z.P. Inelastic analysis of structures. - New York: Wiley, 2002. - 734 p.
14. Bazant Z.P., Zi G. Microplane constitutive model for porous isotropic rocks // Int. J. Numer. Anal. Meth. Geomech. - 2003. - 27:25-47
15. Lai T.Y., Borja R.I., Duvernay B.G., Meehan R.L. Capturing strain localization behind a geosynthetic-reinforced soil wall // Int. J. Numer. Anal. Meth. Geomech. - 2003. - 27:425-451
16. Papamichos E. Constitutive laws for geomaterials, oil and gas science and technology // Rev. IFP. - 1999. -V. 54. - No. 6. - P. 759-771.
17. Borja R.I., Regueiro R.A., Lai T. Y. FE modeling of strain localization in soft rock // J. Geotech. Geoenviron. Engrg., ASCE. - 2000. -V. 126.- P. 335-343.
18. Vardoulakis I. Behavior of granular materials // Handbook ofmaterials behavior models / Ed. by J. Lemaitre. - San Diego: Academic Press, 2001. - V. 3. - P. 1093-1105.
19. Lubarda V.F., Mastilovich S., Knap J. Some comments on plasticity postulates and non-associative flow rules // Int. J. Mech. Sci. - 1996. -V. 38. - No. 3. - P. 247-258.
20. Colmenares L.B., Zoback M.D. A statistical evaluation of intact rock failure criteria constrained by polyaxial test data for five different rocks // Int. J. of Rock Mechanics and Mining Sciences. - 2002. - V. 39. -P. 695-729.
21. Bardet J.P. Lode dependences for isotropic pressure-sensitive elasto-plastic materials // J. of Appl. Mech. - 1990. - V. 57. - P. 498-506.
22. Fairhurst C., Cook N.G.W. The phenomenon of rock splitting parallel to the direction of maximum compression in the neighborhood of a surface // Proc. First Congress International Society for Rocks Mechanics, Lisbon, 1966. - P. 687-692.
23. Germanovich L.N., Salganik R.L., Dyskin A.V., Lee K.K. Mechanisms of brittle fracture of rock with pre-existing cracks in compression // Pageoph. - 1994. - V. 143. - No. 1/2/3. - P. 117-149.
24. Фриденталь А., ГейрингерX. Математические теории неупругой сплошной среды. - М.: Гос. изд. физ-мат. лит., 1962. - 432 с.
25. Черепанов О.И. Численное моделирование деформации материалов с учетом неустойчивой ветви G-e-диаграммы // Физ. мезомех. -1999. - Т. 2. - № 1-2. - С. 5-16.
26. Zhao J. Applicability of Mohr-Coulomb and Hoek-Brown strength criteria to the dynamic strength of brittle rock // Int. J. of Rock Mechanics and Mining Sciences. - 2000. - V. 37. - P. 1115-1121.
27. Huismans R.S., Podladchikov Y.Y., Cloetingh S. Transition from passive to active rifting: Relative importance of astenospheric doming and passive extension of the lithosphere // J. of Geoph. Research. -2001. - V. 106. - No. B6. - P. 11.271-11.291
28. Новожилов В.В. О пластическом разрыхлении // Прикладная математика и механика. - 1965. - Т. 29. - Вып. 4. - С. 681-689.
29. Уилкинс М.Л. Расчет упругопластических течений // Вычислительные методы в гидродинамике. - М.: Мир, 1967. - С. 212-263.
30. Стефанов Ю.П. Локализация деформации и разрушение в геоматериалах. Численное моделирование // Физ. мезомех. - 2002. -Т.5. - №5. - С. 107-118.
31. Ахмадеев Н.Х. Динамическое разрушение твердых тел в волнах напряжений. - Уфа: БФАН СССР, 1988. - 168 с.
32. Майчен Дж., Сак С. Метод расчета «Тензор» // Вычислительные методы в гидродинамике. - М.: Мир, 1967. - С. 185-211.
33. Gulidov A.I., Fomin V.M., Shabalin I.I. Mathematical simulation of fracture in impact problems with formation of fragments // Int. J. of Fracture. - 1999. - V. 100(2). - P. 121-131.
34. Гриднева В.А., Немирович-Данченко М.М. Численный расчет поведения материала образца с трещиной при одноосном растяжении // Механика деформируемого твердого тела. - Томск: ТГУ, 1985. - С. 59-63.
35. Немирович-Данченко М.М. Модель гипоупругой хрупкой среды: применение к расчету деформирования и разрушения горных пород // Физ. мезомех. - 1998. - Т. 1. - № 2. - С. 107-114.
36. Chen Y.M., Wilkins M.L. Stress analysis of crack problems with a three-dimensional time-dependent computer program // Int. J. of Fracture. - 1976. - 12(4). - P. 607-617.
37. Вычислительные методы в механике разрушения / Под. ред. С. Ат-лури. - М.: Мир, 1990. - 392 с.
38. Nishioka T. Computational dynamic fracture mechanics // Int. J. of Fracture. - 1997. - V. 86. - P. 127-159.
39. Bourago N.G. A survey on contact algorithms // Proc. of Workshop «Grid Generation: Theory and Applications» / Ed. by S.A. Ivanenko, V.A. Garanzha. - Moscow: Computing Centre of RAS, 2002. - P. 4259.
40. Бураго Н.Г., Кукуджанов В.Н. Численное решение задач континуального разрушения / Препринт № 746. - М.: Институт проблем механики РАН, 2003. - 38 с.
41. Johnson G.R., Stryk R.A. Symmetric contact and sliding interface algorithms for intense impulsive loading computations // Comput. Meth. Appl. Mech. Eng. - 2001. - 190. - P. 4531-4549.
42. Stefanov Yu.P. Wave dynamics of cracks and multiple contact surface interaction // Theor. and Appl. Frac. Mech. - 2000. - V. 34/2. - P. 101108.
43. 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.
44. Стефанов Ю.П., Поболь И.Л., Князева А.Г., Гордиенко А.И. Рост трещины вблизи границы раздела разнородных материалов в условиях сжатия // Физ. мезомех. - 2002. - №1 . - С. 81-88.
45. Макаров П.В., Трубицын А.А., Трубицына Н.В., Кузнецов П.В., Смолин И.Ю., Стефанов Ю.П., Ворошилов С.П. Экспериментальное и теоретическое исследование разрушения углей и расчет выхода пылевых частиц. II. Численное изучение разрушения угля на мезо- и макроуровнях // Физ. мезомех. - 2004. - Спец. вып. - Ч.2.-С. 249-252.
46. Стефанов Ю.П. Численное исследование поведения упругоидеальнопластических тел, содержащих неподвижную и распространяющуюся трещины, под действием квазистатических и динамических растягивающих нагрузок // Физ. мезомех. - 1998. -№2. - С. 81-93.
47. ПсахьеС.Г, СмолинА.Ю., СтефановЮ.П., МакаровП.В., Чертов М.А. Моделирование поведения сложных сред на основе совместного использования дискретного и континуального подходов // Письма в ЖТФ. - 2004. - Т. 30. - Вып. 17. - С. 7-13.
48. Смолин А.Ю., Стефанов Ю.П., Псахье С.Г. Совместное использование дискретного и континуального методов для моделирования процессов деформации и разрушения в области контактного взаимодействия. // Физ. мезомех. - 2004. - Спец. вып. - Ч. 1 - С. 70-73.
49. Wilkins M.L. Use of artificial viscosity in multidimensional shock wave problems // J. Comput. Phys. - 1980. - V. 36. - P. 281-303.
50. Johnson G.R., BeisselS.R. Damping algorithms and effects for explicit dynamics computations // Int. J. Impact Engineering. - 2001. - V. 25. -P. 911-925
51. Gerbault M., Poliakov A.N.B., Daignieres M. Prediction of faulting from the theories of elasticity and plasticity: What are the limits? // J. Struct. Geolog. - 1998. - V. 20. - No. 2/3. - P. 301-320.
52. Labuz J.F., Dai S.-T, Papamichos E. Plane-strain compression of rock-like materials // Int. J. of Rock Mech. and Min. Sci. and Geomech. Abstr. - 1996. - V. 33. - No. 6. - P. 573-584.
