CC BY
40
УДК 539.3
Применение метода механической аналогии для численного моделирования разрушения керамических композитов ZrO2-Al2O3 в трехмерной постановке
М.О. Еремин
Национальный исследовательский Томский государственный университет, Томск, 634050, Россия Институт физики прочности и материаловедения СО РАН, Томск, 634055, Россия
В работе представлено обобщение метода построения регулярных сеток на трехмерную область на основе уравнений механики деформируемого твердого тела, описывающих деформацию конечного объема. Представлены расчеты разрушения керамических композитов с включениями на основе разработанной модели квазихрупкой среды с учетом накопления повреждений.
Ключевые слова: метод механической аналогии 3D, численное моделирование, разрушение, накопление повреждений, нелинейная динамическая система, режимы с обострением
Numerical simulation of fracture of ZrO2-AhO3 ceramic composites
M.O. Eremin
National Research Tomsk State University, Tomsk, 634050, Russia Institute of Strength Physics and Materials Science SB RAS, Tomsk, 634055, Russia
The regular grid generation method was generalized for a 3D region based on the solid mechanics equations that describe deformation of a finite volume. Test calculations of fracture of ceramic composites with inclusions are represented which are based on the developed model of quasi-brittle media with regard to damage accumulation.
Keywords: 3D mechanical analogy method, numerical simulation, fracture, damage accumulation, nonlinear dynamic system, blowup regimes
1. Введение
В работе [ 1 ] описан метод построения регулярных сеток в двумерной области на основе уравнений, описывающих продольную деформацию пластинки. С использованием данного метода был решен ряд задач моделирования деформации материалов с пористыми покрытиями, например [2]. Метод представляет интерес для задач, где необходимо создавать сетки с учетом сложной геометрии структурных элементов, например границ раздела фаз в механике композиционных материалов. Особенно это важно становится при применении численных схем интегрирования уравнений механики деформируемого твердого тела, использующих шаблоны квадратных (4 узла) ячеек в двумерных задачах или кубических (8 узлов) ячеек в трехмерных задачах. Такие схемы обладают низкой аппроксимационной вязкостью и, как следствие, лучше описывают пространст-
© Еремин М.О., 2015
венную локализацию деформаций в нагружаемых твердых телах, однако зачастую границы раздела фаз в таких схемах представляют собой ломаные линии, что ухудшает гладкость решения. В связи с чем использование системы эллиптических уравнений для получения гладких регулярных сеток на границах раздела фаз представляется целесообразным.
В данной работе метод построения регулярных сеток обобщен на трехмерный случай и получены соответствующие уравнения. Построенные сетки использованы для тестовых расчетов разрушения керамических композитов с одиночным и несколькими сферическими включениями.
2. Описание метода
В отсутствие массовых сил и инерционного члена закон сохранения количества движения выражается тремя уравнениями равновесия:
дахх , дТху + дт хг = 0,
дх ду дг
дт дауу дт
ух 1 уу 1 у- = 0,
дх ду дг
+ да гг = 0.
(1)
Эх ду дz Далее компоненты тензора напряжений а у из обобщенного закона Гука, выраженные через компоненты вектора перемещений и:
а„ =
ауу =
ат =
(1 + а)(1 - 2 а)
Е
(1 + а)(1 - -2 а)
Е
ди (ди
(1 -а)и + а
(1 -а)-
дх
ди.,
у ди. у +
ду
- + а
дz
ди,
Тху
Тх2 =
т у- =
(1 + а)(1 -2 а)
(дих
(1 -а) и + а(
дг
ду дих
дх дг
дих +диу ^
дх ду
диу )
2(1 + а) / ду + дх
Е (дих .д^
2(1 + а) дх
Е (диу ди-+-2-
(2)
)
2(1 + а) / дг ду были подставлены в соответствующие уравнения равновесия. В результате была получена следующая система уравнений, известная как система уравнений Ламе:
д2их 1 - 2а д2их 1 - 2а д2их
- +--X +--X +
дх2
2(1 -а) ду2 2(1-а) дг2
1 д2и у
1 д2и2
= 0,
2(1 -а) дхду 2(1-а) дхдг
д2иу 1 - 2а д2иу 1 - 2а д2иу ду2 + 2(1 -а) дх2 +
1 д2их
2(1 -а) дг2
1 д2и
(3)
2(1 -а) дхду 2(1 -а) дудг
= 0,
д2и2 1-2а д2и2 1 -2а ди
-X +--X +--^ +
дг 2 +
2(1 -а) дх2 2(1 -а) ду2
1 д 2их
1 д2иу
= 0,
2(1 -а) дхдг 2(1 -а) дудг где о — коэффициент Пуассона. Для получения конечно-разностного аналога системы уравнений (2) воспользуемся следующими конечно-разностными аппроксимациями частных производных:
д2 f = А+1,у,к - 2Уг,у,к + /-1,у',к
дх
Ах2
д7 ./г, у+1,к ^./г, у ,к + Уг", у-1,к
ду2 Ау2
д7 ./г, у,к+1 - 2 А, у ,к + Уг, у,к-1
дг 2 Аг2 '
А+1, у+1, к - Уг-1,у+1,к - Уг+1,у-1,к + Уг-1,у-1,к
дхду 4АхАу
д2 А = А+1, у, к+1 - ./г-1,у,к+1 - Уг+1,у,к-1 + ./¿-1,у,к-1
дхдг 4АхАг
д2/ = ./г, у+1,к+1 - Л,у -1,к+1 - Л,у +1,к-1 + ./г,у -1,к-1
дудг 4АуАг
где / ук — некоторая функция, определенная в узлах сетки. Если конечный объем сплошной среды разбит регулярной сеткой с количеством узлов (т + 1)(п + 1)х х(/+ 1), где т, п, I — количество ячеек вдоль осей X, У, Z соответственно, то в результате аппроксимации получаем следующую систему линейных алгебраических уравнений относительно неизвестных перемещений
Цх_г,у ,к' Цу _г,к' _г,у',к :
их_г,у',к = ах (их_г +1,у ,к + Цх_г-1,у',к ) + + Рх (их_г,у +1,к + их_г,у -1,к ) + + Ух (их _ г, у к+1 +их _ г, у ,к-1) + + 8х(иу_г+1,у +1,к-иу_г-1,у +1,к -
- иу _ г+1, у -1,к + иу _ г -1, у -1,к ) +
+ ех (Ц2 _г +1,у ,к+1 - Ц2 _г -1,у,к+1 -
- Ц2_г+1,у,к-1 + Ц2_г-1,у,к-1 ),
1 1 - 2а
(5)
а х =
Ах2 V Рх 2(1 -а)Ау2,
Ух =
1 - 2а
8 х =-
1
е х =-
2(1 -а)Аг2 у 8(1 -а)АхАуЛук
1_
8(1 -а)АхАЛук '
2
Лук = — +
1 - 2а 1 - 2а +-
Ах (1 -а) Ау2 (1 -а) Аг2 '
иу _ г,у,к = ау (иу _ г,у +1,к + иу _ г,у -1,к ) +
+ Ру (иу _ г+1,у ,к + иу _ г-1,у ,к ) +
+ Уу (иу _г,у,к+1 +иу _г,у,к-1) + + 8у (их _ г+1, у +1, к -их _ г-1, у +1, к -
- их_г+1,у-1,к + их_г-1,у-1,к ) +
+ еу (иг_ г,у +1,к+1 -иг_ г,у -1,к+1 -
- иг_г,у+1,к-1 + иг_г,у-1,к-1),
1 Ру ^ 1 - 2а
а у = л 2
А/Вк' '' 2(1 -а)Ах2V
У V =
е V =
1 - 2а
8 V =-
1
2(1-а) Д2 Вт У 8(1 -а)ДхДуВ# 1
% = ТГ+
1-2а 1 - 2а - +-
Д/ (1 - а) Дх2 (1 - а) Д,2 '
_г,у,к = а2 (иг_г,у,к+1 + и,_г,у,к-1 ) +
+ Р, (иг _ г+1, ]Л +и2 _ г -1 м ) +
+ У, (и,_г,у +1,к + _г,у -1,к ) + + (их_г+1,>,¿+1 - их_г-1,>,¿+1 -
- их_г+1,у,к-1 + их_г-1,у,к-1 ) +
+ ег _ г, 7+1,к+1 V _ г, у-1, ¿+1 -
- ^_г,у +1,к-1 + ^V_г,у-1,к-1); 1 - 2а
а, = - "
(7)
' 2(1 -а)Дх2С,
У, =
1 - 2а
8, =-
ук 1
е, = -
2(1 -а)Ду2С# 8(1 - а)ДхД,Сук 1
8(1 -аДДСк'
Су = ТГ+
1 - 2а 1 - 2а - +-
Д^2 (1 -а)Дх2 ' (1 -а)Ду2' Данная система при заданных граничных условиях имеет 3(т - 1)(и - 1)(/ - 1) неизвестных. Как отмечалось в работе [1], система может быть решена при помощи метода Зейделя [3].
3. Математическая постановка задачи разрушения композита с одиночным включением
Для моделирования процесса разрушения керамического композита с одиночным включением будем решать полную систему динамических уравнений механики деформируемого твердого тела, которая включает в себя фундаментальные законы сохранения массы и количества движения, геометрические соотношения (8) и две группы определяющих соотношений (9) и (10)-(12):
Р^ = Р0*0> Р^г = агу,у > Щ = V, у + V у г > Щу = V, у - у; Р = -к(6х -ёР), ёх = е^Г,
¡¡гу + БЛ 6>ку - Бку Щк = 2^|4х - 3 0Х8г,--егу I;
f(ау) = -аР + 4Гг -7, g (ау) = J2 - ЛР(27 +аР) + const, 7 = 7о(1 - В).
(8)
(9)
(10)
Здесь р0, р, У0, V — начальное и текущее значение плотности и объема материала соответственно; V — компоненты вектора скорости; Р — давление, агу, Бу — компоненты девиатора тензора напряжений; еХ — компоненты тензора скорости деформации. Для учета поворота элемента среды как целого используется корота-ционная производная по времени Яуманна, 8гу — символ Кронекера; J2 = 1/2 БуБ у — второй инвариант девиатора тензора напряжений; К — модуль объемного сжатия; Ц — модуль сдвига; а — коэффициент внутреннего трения; Л — коэффициент дилатансии; g (агу) — пластический потенциал, обеспечивающий неассоциированный закон течения (таким образом, процесс дилатансии независим от внутреннего трения); У— сцепление (сдвиговая прочность среды при нулевом давлении), которое уменьшается от начального значения 7 по мере накопления поврежденности D. Применяя основное соотношение теории пластичности е Р = = Лдg (а у у да у , где пластический множитель X определяется из условия удовлетворения напряжений условию текучести (10), получаем следующее соотношение для компонент тензора скоростей пластических деформаций [4]:
а
еу =1 Бу + —Л(7--Jl)8ii X
= еР
77 *
■чу 1 з1 з иг*' ^ (11)
Релаксация напряжений за счет накопления пластических деформаций стабилизирует процесс деформирования и выступает в качестве отрицательной обратной связи.
Функцию деградации среды (поврежденность) D = = D(t, Ца, о) < 1 представим в виде зависимости от инварианта напряженного состояния асиг и вида напряженного состояния Ца (коэффициента Лоде-Надаи):
В = 1
' [Н(Ца)(асиг -а0)2 + (1 -н(ЦаЖа^ -а0)2№
а2[Н(ЦаК + (1 - Н(Ца Ж*]
а* = а0* (1.01 + Ца) , Ца = 2
Б2 Б3 Б1 - Б3
- 1,
(12)
где Н(х) — функция Хевисайда; асиг = аint - аР, аint — интенсивность тензора напряжений; а0, а0 — начальные значения напряжений на упругой стадии, по достижении которых в материале начинается накопление повреждений в областях сжатия и растяжения соответственно, причем а0 << а0. Таким образом, повреждения в областях растяжения-сдвига (Ца < 0) начинают накапливаться при существенно меньших напряжениях, чем при Ца >0 в областях сжатия-сдвига. Скорость накопления повреждений для локальных областей, где Ца < 0, также существенно выше, чем в областях сжатия-сдвига (Ца > 0). Этот процесс дополнительно регулируется параметром а* в (12), формируя в квазихрупкой среде существенно меньшую прочность при растяжении-сдвиге. Таким образом, отклик среды на вид напря-
женного состояния (ее текущая прочность) формируется в среде в процессе ее нагружения. Следовательно, прочностные параметры будут деградировать существенно быстрее в тех областях среды, где коэффициент Лоде-Надаи ца <0, т.е. преимущественный характер напряженно-деф ормированного состояния определяется растяжением-сдвигом. В выражении (12) а 0* — параметр модели; 4 имеет смысл характерного времени процесса разрушения; S1, S2, Sз — главные значения девиатора тензора напряжений. Функция деградации в системе выступает в качестве положительной обратной связи, приводящей к неустойчивости процесса деформирования.
Система (8)-(12) была решена на основе метода Уилкинса, описанного в работе [5]. Прежде чем перейти к результатам моделирования, необходимо несколько слов сказать об описании процесса разрушения, примененном в данной работе. Идеи синергетики (нелинейной динамики) проникли в механику и физику твердого тела. В ряде работ [6-8] деформируемое твердое тело стали рассматривать как синергетическую систему, в которой протекают процессы самоорганизации под приложенными воздействиями. В твердом теле под нагрузкой формируются разномасштабные диссипативные структуры. На микроуровне — это скопления различного рода дефектов, например дислокаций или микротрещин, на мезоуровне — это мезотрещины, мезополосы локализованной пластической деформации, на макроуровне — это магистральные трещины. Таким образом, процесс разрушения является многомасштабным и охватывает всю иерархию масштабов. Переход на следующий масштаб является результатом потери устойчивости на предыдущем. Эта последовательность продолжается до тех пор, пока не будет достигнут макромасштаб и не наступит макроразрушение как глобальная потеря устойчивости. Эта идея является одной из основных
идей физической мезомеханики материалов. Современная механика деформируемого твердого тела объединяет в себе идеи классической механики сплошных сред, физической мезомеханики и нелинейной динамики (синергетики).
Существенный вклад в исследование формирования диссипативных структур в нелинейных динамических системах был в работе [9]. Концепция режимов с обострением является важной и для понимания процесса разрушения твердых тел.
Исходя из идей [9], для корректного описания формирования диссипативных структур в нелинейной динамической системе необходимо описать положительные и отрицательные обратные связи. Эти связи описываются посредством эволюционных уравнений. Для деформируемого твердого тела эволюционные уравнения выражаются системой (8)-( 12). Они объединяют макроуровень — связь скоростей роста/релаксации напряжений со скоростями деформаций, мезоуровень, на котором явно вводятся важнейшие структурные элементы, определяющие ключевые механизмы деформационных процессов в твердых телах, микроуровень, интегрально отражающий вклад микромасштабов в деформационный процесс посредством кинетических соотношений для скоростей генерации дефектов, неупругих деформаций или повреждений [6].
4. Результаты моделирования и обсуждение
В данной работе в модели композита для простоты форма включения была взята сферической (рис. 1). Использована схема преобразования включения в форме куба в сферу, аналогично процедуре преобразования квадрата в круг, примененного в работе [2].
В качестве материалов композиции были использованы диоксид циркония ZrO2 (матрица) и корунд А1203
Рис. 1. Разбиение образца композита с одиночным включением сеткой с применением метода построения регулярных сеток в 3Б (показано центральное сечение модели и увеличенное изображение сетки во включении и вблизи него)
Таблица 1
Физико-механические свойства композита, использованные при моделировании
Материал Плотность р0, г/см3 Модуль объемного сжатия К, ГПа Модуль сдвига т, ГПа Сцепление У0, МПа Коэффициент внутреннего трения а Коэффициент дилатансии Л
ZrO2 5.700 143 66.15 2100 0.50 0.10
А12°3 3.984 346 160.00 3740 0.55 0.12
(включение). Физико-механические свойства материалов композиции представлены в табл. 1. При моделировании были использованы следующие параметры разностной сетки: количество узлов по осям координат Nx = 200, NV = 200, Nz = 200, объем модели композита составляет 50x50x50 мкм3. Нагружение осуществлялось приложением постоянной скорости смещения узлов сетки вдоль оси Oz на нижней и верхней гранях образца, ортогональных оси Oz.
При одноосном сжатии композита структурная неоднородность приводит к формированию неоднородного напряженно-деформированного состояния вокруг включения (рис. 2). На рис. 2 представлено распределение инварианта напряженного состояния: асиг = = а^ - аР. Поскольку предельная поверхность, ограничивающая упругую область напряженного состояния, записывается в форме Мизеса-Шлейхера, представляющей конус в пространстве напряжений, то анализировать напряженно-деформированное состояние удобно именно в такой форме. Области концентрации напряжений отмечены окружностями на рис. 2. Именно в этих
областях начинают зарождаться полосы локализованных повреждений в матрице, что подтверждается распределением инварианта напряженного состояния, соответствующего следующему моменту времени. Происходит существенная релаксация напряжений в матрице на первоначальных концентраторах, уровень напряжений уменьшается в 5 раз, в то же время в вершине растущих поверхностей локализованных повреждений наблюдается концентрация напряжений, почти в 1.52.0 раза превышающих уровень напряжений на начальных концентраторах, обеспечивающая возможность дальнейшего роста поверхностей.
Поля, отражающие вид напряженного состояния по параметру Лоде-Надаи, также показывают кардинальное его изменение — от сжатия со сдвигами (до стадии начала формирования областей локализованных повреждений (Ца > 0)) до растяжения со сдвигами для следующего момента времени (Ца < 0), когда на интерфейсе зарождаются области локализованных повреждений.
По мере накопления повреждений в среде еще на упругой стадии, т.е. повреждений меньших масштабов,
0 Концентраторы напряжений
ш
в
*
а, МПа -490
-400
300
'-'-228
а, МПа -710
-600
-400
200
1-1-34
Ца 1-0-800
-0.400
-0.000
-0.400 -0.673
Рис. 2. Распределение инварианта напряженного состояния и вида напряженного состояния по параметру Лоде-Надаи для двух последовательных времен г1 (а) и г2 (б) в осевом сечении образца
I—0.00
Рис. 3. Эволюция поврежденности композита в процессе деформирования
чем мезо и макро, предельная поверхность, ограничивающая упругую область, постепенно сужается, и скорость сужения будет постоянно увеличиваться за счет действия положительной обратной связи. Когда значение инварианта напряженного состояния асиг сравняется с текущей предельной поверхностью, в материале начнется стадия активного разрушения, или процесс разрушения перейдет на мезомасштабный уровень. Когда в материале сформируются магистральные полосы локализованной неупругой деформации, это будет означать, что процесс разрушения достиг макромас-штабного уровня. За счет действия положительной обратной связи на заключительных этапах деформирования процесс разрушения развивается в режиме с обострением, и происходит глобальная потеря устойчивости всей системы.
На рис. 3 показаны этапы процесса разрушения композита с одним включением. Расчеты показывают, что при осевом сжатии композита в первую очередь происходит отслоение матрицы от включений. Объясняется это наличием растягивающих напряжений на интерфейсах, расположенных параллельно оси сжатия (на рис. 3 ось сжатия ориентирована вертикально). Кроме системы формирующихся субвертикальных поверхностей локализованных повреждений, расположенных в основном вблизи центральной оси образца, параллельной оси сжатия, в образце также формируются поверхности локализованных повреждений, ориентированные вдоль максимальных касательных напряжений.
Характер разрушения композита является хрупким, на что указывает не только распределение поврежден-ности в объеме, но также и вид макроскопической кривой нагружения (рис. 4), на которой отчетливо виден незначительный этап неупругого деформирования, выражающийся отклонением диаграммы от упругого закона. На заключительном этапе деформирования процесс разрушения развивается в режиме с обострением, на диаграмме этот этап представлен ниспадающей веткой. Данная кривая нагружения получена при следующих параметрах модели накопления повреждений: а0 = = 0.270,а а0 = 0.0570, а0* = 0.1, £ = 10-9, tZ = 10
к-7
Рис. 4. Макроскопическая кривая нагружения композита в координатах «осевое напряжение - осевая деформация»
0
•4 г 2
0
2 X
0 0
в объеме композита
На рис. 5 показана динамика формирования поверхностей локализованных повреждений в объеме композита, содержащего порядка 20 сферических включений разных размеров.
5. Заключительные замечания
В результате обобщения метода построения регулярных сеток на случай 3D получена система линейных алгебраических уравнений относительно неизвестных перемещений узлов сетки при заданных граничных условиях. Метод применен для построения моделей композита с одиночным и несколькими включениями,
имеющими сферическую форму. В силу малости предполагаемых деформаций в системе уравнений Ламе метод имеет ограничения в области построения сеток с достаточно сложной геометрией границ через единое преобразование всей области. Однако для построения относительно гладких границ раздела метод применим, что подтверждают представленные в статье результаты расчетов.
Полученные ранее результаты о возникновении растягивающих напряжений на интерфейсах [2, 10] в данной работе подтверждаются результатами моделирования осевого сжатия композитов с гладкими границами раздела в трехмерной постановке. На интерфейсах образуются области локализованных повреждений, где меняется вид напряженного состояния от сжатий со сдвигами на начальных этапах деформирования до растяжений со сдвигами на более поздних этапах. Дальнейшее деформирование композита приводит к формированию развитой сети поверхностей локализованных повреждений, которая достигает макроуровня на ниспадающей ветви макроскопической кривой нагружения.
Для корректного описания особенностей хрупкого/ квазихрупкого разрушения и перехода процесса разрушения к сверхбыстрому этапу эволюции — режиму с обострением — в системе уравнений прописываются положительные и отрицательные обратные связи [11]. Полная система уравнений механики деформируемого твердого тела с учетом кинетических уравнений для накопления повреждений позволяет качественно и количественно верно описать процесс разрушения композитов на различных масштабных уровнях.
Характер разрушения оказывается самоподобным, начиная с масштабов элементарной ячейки и заканчивая масштабом всего образца в целом. Это подтверждается интегральной оценкой характера накопления повреждений на различных масштабных уровнях. Накопление повреждений в отдельно взятой ячейке, мезообъеме или образце в целом развивается схожим образом — наблюдается медленная квазистационарная стадия и затем наступает режим с обострением. На микро- и мезоуров-нях переход в режим с обострением на соответствующем уровне сопровождается локальной потерей устойчивости в элементарной ячейке и мезообъеме соответственно. На макроскопическом уровне переход в режим с обострением процесса разрушения сопровождается глобальной потерей устойчивости всего образца и обрывом диаграммы деформирования. Подобное поведение свойственно реальным нагружаемым твердым телам, в которых процесс разрушения охватывает всю иерархию масштабных уровней и является самоподобным [12].
Работа выполнена при поддержке программы фундаментальных исследований СО РАН № Ш.23.2.3 и в рамках программы повышения конкурентоспособности ТГУ.
оУ-— л 1 1 1 1 Г
Литература
1. Игнатьев А.А. Построение регулярных сеток с помощью механической аналогии // Математическое моделирование. - 2000. -Т. 12. - № 2. - С. 101-105.
2. Balokhonov R., Zinoviev A., Romanova V., Martynov S. A numerical simulation of the deformation and fracture of a material with a porous polysilazane coating // Int. Conf. on Physical Mesomechanics of Multilevel Systems 2014: AIP Conf. Proc. - 2014. - V. 1623. - P. 51-54 (doi: 10.1063/1.4898880).
3. Ханова А.А. Численное решение уравнений и систем уравнений /
Под ред. В.В. Лаптева. - Астрахань: Изд-во АГТУ, 2001. - 27 с.
4. Stefanov Yu.P., Chertov M.A., Aidagulov G.R., Myasnikov A.V. Dynamics of inelastic deformation of porous rocks and formation of localized compaction zones studied by numerical modeling // J. Mech. Phys. Solids. - 2011. - V. 59. - P. 2323-2340.
5. Wilkins M.L. Computer Simulation of Dynamic Phenomena. - BerlinHeidelberg: Springer-Verlag, 1999. - 246 р.
6. Макаров П.В. Эволюционная природа блочной организации геоматериалов и геосред. Универсальный критерий фрактальной делимости // Геология и геофизика. - 2007. - Т. 48. - №2 7. - С. 724746.
7. Панин В.Е. Синергетические принципы физической мезомеханики
// Физ. мезомех. - 2000. - Т. 3. - № 6. - С. 5-36.
8. Пантелеев И.А., Плехов О.А., Наймарк О.Б. Нелинейная динамика
структур обострения в ансамблях дефектов как механизм формирования очагов // Физика Земли. - 2012. - № 6. - С. 43-54.
9. Курдюмов С.П. Режимы с обострением. Эволюция идеи. - М.: Физматлит, 2006. - 312 с.
10. Romanova V.A., Balokhonov R.R., Schmauder S. The influence of the reinforcing particle shape and interface strength on the fracture behavior of a metal matrix composite // Acta Mater. - 2009. - V. 57. -P. 97-107.
11. Смолин И.Ю., Еремин М.О., Макаров П.В., Буякова С.П., Кульков С.Н. Численное моделирование механического поведения модельных хрупких пористых материалов на мезоуровне // Вестник ТГУ. Математика и механика. - 2013. - № 5. - С. 78-90.
12. Пантелеев И.А., Froustey C., Наймарк О.Б. Структурно-скейлин-говые переходы и автомодельные закономерности развития землетрясений // Вычислительная механика сплошных сред. - 2009. -Т. 2. - № 3. - С. 70-81.
Поступила в редакцию 21.01.2015 г., после переработки 07.04.2015 г.
Сведения об авторе
Еремин Михаил Олегович, к.ф.-м.н., мнс ТГУ, мнс ИФПМ СО РАН, eremin@ispms.tsc.ru
