CC BY
91
УДК 539.37
Моделирование пластической деформации на основе теории ортотропного континуума Коссера
В.М. Садовский1'2, М.А. Гузев34, О.В. Садовская1, Ch. Qi4
1 Институт вычислительного моделирования СО РАН, Красноярск, 660036, Россия 2 Сибирский федеральный университет, Красноярск, 660041, Россия 3 Институт прикладной математики ДВО РАН, Владивосток, 690041, Россия 4 Beijing University of Civil Engineering and Architecture, Beijing, 100044, China
К анализу пластической деформации структурно-неоднородных материалов применяется метод прямого численного моделирования на основе теории ортотропного упругопластического континуума Коссера с условием пластичности, учитывающим как сдвиговый, так и ротационный характер необратимой деформации. В рамках предположения о блочной структуре материала с упругими блоками, взаимодействующими через податливые пластические прослойки, это условие ограничивает касательные компоненты несимметричного тензора напряжений, характеризующие сдвиги, а также моментные напряжения, достижение которыми предельных значений приводит к необратимому изменению характеристик кривизны деформированного состояния континуума. Уравнения поступательного и вращательного движения совместно с определяющими соотношениями модели формулируются в виде вариационного неравенства, корректно описывающего как состояние упругопластического деформирования при активном нагружении, так и состояние упругой разгрузки. Для численной реализации математической модели используются параллельный вычислительный алгоритм и авторский программный комплекс для многопроцессорных вычислительных систем кластерной архитектуры. Исследуется задача о сдавливании прямоугольного массива блочной горной породы типа кирпичной кладки шероховатой недеформируемой плитой, совершающей равноускоренное вращательное движение. Рассмотрено влияние пределов текучести податливых прослоек при сдвиге и при изгибе на напряженно-деформированное состояние массива. Наряду с полями перемещений, напряжений, моментных напряжений и угла поворота структурных элементов анализируется поле пластической диссипации энергии в массиве. Полученные результаты могут служить обоснованием гипотезы о преимущественном влиянии кривизны на образование зон локализации пластической деформации на мезоуровне в материалах с микроструктурой.
Ключевые слова: микроструктура, моментные напряжения, континуум Коссера, упругость, пластичность, вариационное неравенство, высокопроизводительные вычисления
DOI 10.24411/1683-805X-2019-12005
Modeling of plastic deformation based on the theory of an orthotropic Cosserat continuum
V.M. Sadovskii12, M.A. Guzev34, O.V. Sadovskaya1, and Ch. Qi4
1 Institute of Computational Modeling SB RAS, Krasnoyarsk, 660036, Russia 2 Siberian Federal University, Krasnoyarsk, 660041, Russia 3 Institute of Applied Mathematics FEB RAS, Vladivostok, 690041, Russia 4 Beijing University of Civil Engineering and Architecture, Beijing, 100044, China
In the paper, the plastic deformation of heterogeneous materials is analyzed using a direct numerical simulation method based on the theory of an elastic-plastic orthotropic Cosserat continuum, with the plasticity condition taking into account both the shear and rotational mode of irreversible deformation. With the assumption of a block structure of a material with elastic blocks interacting through pliable plastic interlayers, this condition imposes constraints on the tangential components of the asymmetric stress tensor, which characterize shearing, and on the couple stresses, which irreversibly change the curvature characteristics of the deformable state of the continuum as they attain critical values. The equations of translational and rotational motion together with the constitutive equations of the model are formulated as a variational inequality, which correctly describes both the state of elastic-plastic deformation under applied loading and the state of elastic unloading. The numerical implementation of the mathematical model is performed using a parallel computation algorithm and an original software package for cluster-based multiprocessor systems. The developed approach is applied to solve the problem of compressing a rectangular blocky rock mass with brickwork setting by a rough nondeformable plate that rotates with constant acceleration. We study how the yield strength of the pliable interlayers influences the stress-strain state of the rock mass in shear and bending. The field of plastic energy dissipation in the rock mass is analyzed along with the fields of displacements, stresses, couple stresses, and angle of rotation of structural elements. The obtained results can be used to validate the hypothesis about the predominant effect of curvature on plastic strain localization at the mesolevel in materials with microstructure.
Keywords: microstructure, couple stresses, Cosserat continuum, elasticity, plasticity, variational inequality, high performance computing
© Садовский В.М., Гузев М.А., Садовская О.В., Qi Ch., 2019
1. Введение
Начиная с работ Треска и Надаи, в экспериментах по квазистатическому нагружению образцов было многократно показано, что зоны локализации деформаций (линии скольжения, или линии Чернова-Людерса) в пластичных материалах, не обладающих выраженной кристаллической структурой, ориентированы по площадкам максимальных касательных напряжений, и что процесс локализации обусловлен необратимыми сдвигами вдоль этих площадок. Современная ситуация характеризуется осознанием необходимости установления «мостов» между различными подходами в описании процессов пластической деформации на разных структурных уровнях. В связи с этим следует привести работы [1-3], в которых показано, что при изучении механизмов пластичности, определяемых деформационными дефектами, существенное значение для понимания структурных изменений в зоне локализации имеет учет ротационных характеристик кристаллов. Неоднородность на уровне кристаллической микроструктуры служит причиной расхождения с описанием классической теории, в которой не учитываются вращательные степени свободы. В работе [4] развита новая физическая парадигма, учитывающая вихревой характер пластической деформации в твердых телах.
Вращательные степени свободы проявляют себя на уровне блочной структуры в геоматериалах с пластическими свойствами [5, 6]. В них при образовании зон локализации под действием статических нагрузок, превышающих предел упругости, существенную роль играют независимые вращения блоков и приобретаемая за счет градиентов поля вращений кривизна деформированного состояния. Другим примером являются хрупкие горные породы, не обладающие пластичностью, в которых в условиях сильного бокового сжатия формируются зоны локализации деформаций типа домино-структур из вращающихся блоков, запускающие механизм роста трещин сдвига со сверхсейсмической скоростью [7, 8]. Таким образом, экспериментальные исследования указывают на необходимость учета вращательных степеней свободы при математическом моделировании локализации деформаций в структурно-неоднородных материалах.
Описание деформационных процессов с помощью модели континуума Коссера впервые было предложено более 100 лет назад [9]. Было установлено [10], что эта модель, сохраняя связь с теорией дислокаций, обладает новыми свойствами, поэтому представляется перспективным ее применение при моделировании пластической деформации материалов с учетом их микроструктуры. В данной работе анализ локализации пластической деформации выполнен методом численного моделирования в рамках теории континуума Коссера с условием пластичности для напряжений и моментных напряжений.
2. Уравнения континуума Коссера
Простой способ моделирования микроструктуры материала дает приближение блочной среды, образованной прямоугольными упругими блоками одинакового размера с тонкими податливыми межблочными прослойками. В результате применения процедуры осреднения в предположении о малости размера блоков на основе такого представления получается модель орто-тропного континуума Коссера. В зависимости от свойств прослоек континуум может быть упругим, вяз-коупругим или упругопластическим.
Уравнения упругого континуума Коссера в случае плоской деформации имеют следующий вид [11]:
РоV1 = + а12,2' РоV2 = СТ21,1 + а22,2'
J0(Oз = Ц31Д + ^32,2 + Ст21 - ^12 ,
СТ11 = a1V1,1 + bv2,2' СТ 22 = a1V2,2 + ^1,1'
Ст21 = a2(v2,1 -(3) + b2(v1,2 + ®зХ (1)
СТ12 = a2 (v1,2 + (3 ) + b2 (v2,1 - (3 X
й31 = а2(3,1' й32 = а2®3,2-
Уравнения записаны в декартовых координатах относительно проекций вектора линейной скорости v j, отличной от нуля проекции вектора угловой скорости ю3, компонент тензора напряжений Стд и моментных напряжений йд (j, k = 1,2). Здесь р0 — плотность континуума; J0 — произведение момента инерции частицы (блока) на число частиц в единице объема; a1, a2, b1, b2 и a2 — модули упругости ортотропного материала с плоскостями симметрии, параллельными координатным плоскостям. Упругие свойства континуума считаются одинаковыми в направлении обеих осей координат, что соответствует модели блочной среды с квадратными блоками. Точкой над символом обозначается производная по времени, нижним индексом после запятой — производная по соответствующей пространственной переменной x1 или x2.
Определяющие уравнения для скоростей напряжений могут быть представлены в виде:
^СТ11 -^СТ22 = v1,1, ^СТ22 -^СТц = v2,2,
a2CT21 -J^n = V2,1 - (3 > ОгСТСТ 12 - Ь2СТ21 = v1,2 + (3 >
где a1, a2, b1, b2 — модули упругой податливости материала:
a, =
-b, =-
1 2 и2'"1 2
a1 - b a1 -1\
a2 =
2 2 j 2 , b2 2 7 2'
a2 - b2 a2 - b2
Это позволяет представить уравнения (1) в матричной форме:
л ди=+В2 ди+QU
д t
д x1
(2)
с симметричными матрицами-коэффициентами А, 1 2
В, В и антисимметричной матрицей Q. При выполнении неравенств
а > IЬ1 а2 > | Ьг|, «2 > 0
система уравнений (2) относится к классу симметрических ¿-гиперболических систем [12] и систем термодинамически согласованных законов сохранения [13,
14].
Система (2) формулируется относительно вектор-функции и с компонентами V, г2, ю3, стп, ст22, ст21,
ст12, ц31, ц32, записанными в столбец. Характеристи-
1 2
ческая матрица системы Л = п1 В +п2В -сА имеет следующий вид:
-Р0 с 0 0 п1 0 0 п2 0 0
0 -Р0 с 0 0 п2 п1 0 0 0
0 0 - J0 с 0 0 0 0 п1 п2
п1 0 0 -а1 с Ь1 с 0 0 0 0
0 п2 0 Ь1 с -а1 с 0 0 0 0
0 п1 0 0 0 -а2 с Ь2 с 0 0
п2 0 0 0 0 Ь2 с -а2 с 0 0
0 0 п1 0 0 0 0 -с/ а2 0
0 0 п2 0 0 0 0 0 -с/ а.
щую математическую модель упругопластического континуума Коссера. Такая модель формулируется в виде вариационного неравенства
(
(и - и) •
а ди - В ди - в2 ди - Qи
дt Эх, Эх2
Л
> 0,
(3)
Эх1 ил2 и, и е F.
Здесь F — множество допустимых вариаций искомого вектора и; йй — произвольный элемент F. Вариационное неравенство представляет собой формулировку принципа Мизеса максимальной мощности пластической диссипации. Множество F не содержит ограничений на вектор скорости и угловую скорость, граница
Подбирая в ней определенным образом параметры с, п1 и п2, можно выразить симметричные матрицы-коэффициенты системы (2). Например, матрица А получается при с = -1, п1 = п2 = 0, матрица В1 — при п1 = 1, с = п2 = 0. Антисимметричная матрица Q = (qjk) целиком заполнена нулями, кроме ее коэффициентов д36 =
= -Чб3 =1 и ?37 =- ?73 = -1-
Характеристическое уравнение detЛ = 0 имеет шесть нетривиальных корней с = ±с1, с = ±с2 и с = ±с3, которые определяют три скорости слабых ударных волн, распространяющихся в направлении единичного вектора (п1, п2):
V,
Р0
^12 = а1 + а, ±
1
(п2 -п,2)2 + (Ь + ь,)2п2п22,
с1, с2, с3 — скорости продольных, поперечных волн и волн вращательного движения; А12 — трехкратный нулевой корень, соответствующий контактным разрывам.
3. Модель упругопластического континуума
Применяя подход [15, 16], можно на основе системы уравнений теории упругости (2) построить более об-
этого множества в пространстве напряжений и момент-ных напряжений является поверхностью текучести материала. Если вектор и лежит внутри F, то в силу произвольности вариации из (3) следует система уравнений (2), описывающая упругий процесс. Если и — граничная точка, то выполняется ассоциированный закон течения, в соответствии с которым вектор пластической скорости деформации направлен по внешней нормали к границе F.
Важнейшая особенность формулировки модели в виде вариационного неравенства (3) заключается в том, что в ней автоматически учитывается условие упругой разгрузки. Если в некоторой точке континуума после необратимой деформации происходит разгрузка, т.е. если вектор и в фазовом пространстве перемещается с границы множества F вовнутрь множества, то в этой точке выполняется система уравнений теории упругости (2).
Так как в соответствии с принятой моделью поведение ортотропного континуума практически полностью определяется деформационными свойствами податливых прослоек блочной структуры, будем в дальнейшем при моделировании кусковатых горных пород использовать обобщенное условие пластичности типа Кулона-Мора:
1 ^31| ^0 -кцСТ11' 1 ^32| ^0 -Кц°22-
В этом условии т0 и ц0 — пределы текучести материала прослоек при сдвиге и при изгибе; кт и к^ — феноменологические коэффициенты, учитывающие повышение упругих свойств под действием сжимающих напряжений. Нормальные напряжения стп и ст22 в (4) участвуют в качестве неварьируемых параметров упрочнения материала. При наличии положительных (растягивающих) нормальных напряжений условие (4) имеет смысл только до момента отрыва, когда пределы текучести, убывая с ростом напряжений, становятся равными нулю.
Переход в пластическое состояние происходит по мере достижения одного или нескольких равенств в (4). На рис. 1 приведены две схемы, поясняющие смысл условия (4). На первой из них показано, как за счет вращения частиц в прослойках возникают силы реакции, нарушающие симметрию тензора напряжений. Согласно этой схеме, касательные напряжения в континууме Коссера обусловлены двумя факторами — деформацией сдвига и поворотом частиц. Оба эти фактора влияют на процесс перехода в пластичность, поэтому в (4) независимо присутствуют два касательных напряжения. Если поле поворотов неоднородно, то, как показано на второй схеме, за счет линейно распределенной деформации в прослойках возникают моментные напряжения изгиба, которые также учитываются в (4). Достижение моментным напряжением предела текучести приводит к необратимой изгибной деформации про-
0
слойки и необратимому относительному повороту разделяемых ею блоков.
4. Алгоритм корректировки решения
Численное решение краевых задач для вариационного неравенства (3) осуществлялось с помощью процедуры расщепления по физическим процессам и по пространственным переменным. Расщепление по физическим процессам получается автоматически как результат аппроксимации производной по времени:
(
(и - и)
Аи-и - ^ ЭЦ - в 2 эц - QU Дt Эх1 дх2
Л
> 0,
Рис. 1. Касательные напряжения, обусловленные вращением блоков (а), схема моментных взаимодействий (б)
и, и е F.
Здесь и — искомое решение на нижнем временном слое; и — решение на верхнем слое; и — произвольная допустимая вариация О. Решение этого неравенства строится на каждом шаге по времени в два этапа. На первом этапе решается краевая задача для системы дискретных уравнений
а^ = *1 ЭЦ+д 2 ди+QU,
Д t Эх1 дх2 которая аппроксимирует систему уравнений динамической теории упругости (2). На втором этапе полученный при ее решении вектор и корректируется в каждой ячейке пространственной сетки в соответствии с вариационным неравенством
(и - 0) • А (и-и) > 0,Ц и е F. (5)
На этом завершается второй этап.
В силу (5) приращение энергии пластической диссипации за время Дt вычисляется по формуле
д d = и • А (и - и).
Эта формула в точности отвечает интерпретации вариационного неравенства (5) как дискретного варианта принципа максимальной мощности диссипации Ми-зеса.
Так как матрица А симметрична и положительно определена, то решением вариационного неравенства (5) служит проекция вектора и на множество F по евклидовой норме ||и||А = Vи • Аи. Учитывая структуру ограничений (4), определяющих гиперкуб в пространстве напряжений ст21, ст12, ц31, ц32, для нахождения проекции применялся метод простой итерации, основанный на выделении диагонали матрицы А. Фактически в этом алгоритме норма, по которой определяется проекция, заменяется нормой ||и||р = уЦ-Ри, где D— диагональная матрица, после чего проекция строится в явном виде.
В терминах вариационных неравенств итерационный алгоритм записывается так (п — номер итерации):
(и - ип+1) • (рип+1 + (а - р)Цп - АЦ) > 0,
и, и е F.
Вопросы сходимости итераций можно обосновать, используя математический аппарат, изложенный, например, в [16].
На этапе решения упругой задачи применялся вариант двуциклического расщепления по пространственным переменным, в котором на каждом шаге по времени (г, г+Дг) решается следующая серия одномерных задач:
А = В- Ц; и <!■(,) = и"»(,),
дt дх1
, эи(2) _2 Эи(2)
dt
.BU(3)
- = B2
Эх2
U (2)(t ) = U (j)(t + At/ 2),
= QU(3), U (3)(t ) = U (2)(t + At/ 2), (6)
, U (4)(t + At/ 2) = U (3)(t + At ),
t
.BU(4) _2 dU(4)
dt
- = B2
Эх2
АЭи(5) = В1 —, и(5) (t + Д/ 2) = и(4) (t + Д).
Эt Эх1
Вектор-функция и(0) определяется из расчета предыдущего шага по времени, а при г = 0 — из начальных данных задачи. Решение на новом временном слое находится по формуле и ^ + Дt) = и+ Дt).
Краевые задачи для одномерных систем (6) с заданными граничными условиями решаются на основе разностной схемы предиктор-корректор, построенной по принципу схем Годунова. Полученная таким образом схема расщепления устойчива по отношению к ошибкам округления при выполнении одномерного условия устойчивости Куранта-Фридрихса-Леви: с1 Д^ 2 < Дх, где Дх — шаг по пространственной переменной, с1 = у] а1 / р0 — максимальная скорость продольных упругих волн. Схема сохраняет второй порядок точности, если для решения одномерных задач применяется схема второго порядка.
Системы обыкновенных дифференциальных уравнений на третьем этапе расщепления решаются на основе разностной схемы Кранка-Николсон, которая обладает вторым порядком точности, безусловно устойчива и сохраняет механическую энергию системы. Эти свойства играют исключительно важную роль, поскольку уравнения системы содержат малый параметр — объемную плотность момента инерции блочной среды J0,
из-за которого проблема численного исследования модели становится нетривиальной.
Параллельная версия алгоритма для решения вариационного неравенства (3) на многопроцессорных вычислительных системах кластерной архитектуры реализована в виде компьютерной программы, код которой написан на алгоритмическом языке Fortran с использованием основных процедур библиотеки MPI. Детальное описание алгоритма и структуры компьютерной программы приведено в [11, 16, 17].
5. Результаты расчетов
В расчетах использовались приведенные в таблице параметры континуума Коссера для материала типа кирпичной кладки с упругими блоками 0.1x0.1 м, с податливыми межблочными швами заданной толщины 8 [11]. Пределы текучести швов при сдвиге т0 и при изгибе Ц0 были приняты равными 0.86 МПа и 8.8кПа-м соответственно. Упрочнение материала за счет сжатия не учитывалось. Эти значения на два порядка ниже максимального уровня напряжений при расчете упругой задачи. В действительности, соотношение пределов текучести, которое, в первую очередь, зависит от механических свойств материала швов, может изменяться в широком диапазоне, поэтому чтобы проанализировать влияние на процесс деформирования различных факторов (деформаций сдвига и кривизны), одному из двух пределов текучести присваивалось значение, превышающее уровень напряжений в упругой задаче.
Расчеты проводились на вычислительных кластерах ИВМ СО РАН (Красноярск) и МСЦ РАН (Москва). Использовалась достаточно мелкая сетка из 500 x 250 ячеек, размер которых согласован с характерным линейным масштабом микроструктуры материала r0 = V6J0/Р0 ~0.1 м, в том смысле, что на этот масштаб приходится примерно пять ячеек сетки. Для достоверности расчеты, результаты которых приведены на рис. 5 и 6, повторяли на вдвое более мелкой сетке. Принципиальные изменения не были обнаружены.
Рассматривалось монотонное по времени жесткое нагружение прямоугольного массива среды высотой h = 5 м и длиной l = 2h, с линейным профилем скорости на верхней границе v2 = -е0 xj t. Такой профиль соответствует сдавливанию массива недеформируемой тя-
Таблица
Механические параметры континуума Коссера
8, мм р0, кг/м3 J0, кг/м aj, ГПа ôj, ГПа a2, ГПа b2, ГПа а2, МН
0.1 3690 6.15 44.5 19.1 14.8 6.77 22.20
1.0 3650 6.05 38.6 16.5 6.46 4.05 9.60
5.0 3470 5.59 23.8 10.2 2.33 0.53 3.39
Рис. 2. Схема нагружения
желой плитой, совершающей равноускоренное вращательное движение с центром вращения в левом верхнем углу (рис. 2). Выбор программы нагружения основывался на том, что при неравномерном (близком к линейному) распределении давления на массив в нем развивается вихревой механизм пластичности.
Время действия нагрузки 10 = 0.05 с и угловое ускорение плиты е0 = 40 с-2 выбирались для условий квазистатического деформирования, чтобы на решение задачи в конечный момент времени t0 слабо влияли динамические эффекты, связанные с распространением волн. Остальные граничные условия на верхней границе массива формулировались в терминах скоростей:
= 0, ю3 = -е01. На нижней границе ставились кинематические условия жесткой заделки: ,и1 = с2 =ю3 = 0. Левая и правая границы считались свободными поверхностями с нулевыми граничными напряжениями
СТП, СТ21 и Ц31.
Предельное решение статической задачи с такими граничными условиями при высокой сдвиговой пластичности, когда т0 ^ 0, а ц0 превышает уровень упругих моментных напряжений, описывается уравнениями
СТ21 = СТ12 = 0, СТ11,1 = ^22,2 = 0, М-31,1 + ц32,2 = 0,
= = (7)
ц31 = а2Ф3,1, ц32 = а2Ф3,2, где Ф3 — угол поворота частиц. Следовательно, ст22 зависит только от х1 и по физическому смыслу задачи представляет собой линейную функцию ст22 = -ух1, а11 = 0, а ф3 является гармонической функцией, удовлетворяющей двумерному уравнению Лапласа ф311 + + ф3 22 = 0 и смешанным граничным условиям: ф3 = 0 на нижней границе, ф3 = -е012/2 на верхней границе и ф31 = 0 на боковых поверхностях. Соответствующее решение уравнения Лапласа ф3 = - х2 е012/(2И) описывает линейное распределение угла по высоте массива. Интегрирование определяющих уравнений
и1 ,1 = ^^ и2,2 = «2a22, и2,1 -Ф3 = и1 ,2 + Ф3 = 0 с учетом граничных условий позволяет определить коэффициент у и проекции вектора перемещений:
2
Перемещения определяются с точностью до движения тела как жесткого целого и тем не менее в этом решении невозможно удовлетворить граничному условию для горизонтального перемещения на верхней и на нижней границах массива, поскольку при контакте с шероховатой плитой и с основанием возникают локализованные пограничные слои.
Другое предельное решение статической задачи при высокой пластичности, связанной с необратимым изменением характеристик кривизны деформированного состояния, когда ц0 ^ 0, а т0 превышает уровень упругих напряжений, подчиняется уравнениям
Й31 =Й32 = 0 а21 =а12. (8)
В этом случае для перемещений, нормальных и касательных напряжений выполняется система уравнений плоской задачи ортотропной теории упругости с симметричным тензором напряжений. Из системы следует обычное уравнение для угла поворота: Ф3 = = (и21 - и12 )/2, показывающее в силу граничных условий на верхней и нижней границах, что поле поворотов не может быть постоянным во всем объеме среды. Кроме того, отсюда следует существование пограничных слоев в зонах контакта с плитой и основанием. Производные Ф31 и Ф3 2, являются компонентами тензора
-0.040
-0.027
-0.013
0.000
0.013
0.027
0.040
У
и1 =— (Ь1х1 + а1х2), и2 =-уа1х1х2, у = —0т
2а^
Рис. 3. Линии уровня горизонтального (а), вертикального перемещения (б) и угла поворота (в)
кривизны в рассматриваемой системе координат [18]. При учете пластичности они представляются в виде суммы упругих составляющих, пропорциональных мо-ментным напряжениям, и пластических составляющих. В предельном случае, когда ц31, ц32 ^ 0, упругие составляющие исчезают, поэтому характеристики кривизны полностью определяются пластическими процессами.
Численное решение задачи с отличными от нуля, но достаточно низкими пределами текучести повторяет качественные особенности решений уравнений (7) и (8). Это касается как кинематических характеристик (вектора перемещений и угла поворота), так и статических (напряжений и моментных напряжений). Линии уровня горизонтального перемещения на рис. 3, а в первом приближении являются эллипсами, а вертикального перемещения на рис. 3,6 — гиперболами (здесь и далее по горизонтали указана координата х1, а по вертикали— х2). Распределение угла поворота по высоте массива оказывается близким к линейному распределению (рис. 3, в) в основной части за исключением узких зон вблизи верхней и нижней границ.
На рис. 4 представлены линии уровня энергии пластической диссипации dт для сдвиговой пластичности с низким пределом текучести т0 и высоким пределом ц0. Сопоставляются среды с относительно тонкими (8 = 0.1 мм) и толстыми (8 = 5 мм) швами. В случае тонких швов энергия пластической диссипации распределяется практически как в однородном материале. С увеличением толщины происходит перераспределение, в котором проявляется структурная неоднородность материала.
Дополнительные расчеты показали, что при понижении значения ц0 на 2-3 порядка над уровнем упругих
кДж/м 1.00 0.83 0.67 0.50 0.33 0.17 0.00
Рис. 5. Линии уровня энергии пластической диссипации за счет необратимого изменения тензора кривизны
напряжений картина распределения практически не меняется. При этом вклад моментных напряжений dц в полную диссипативную энергию d = dт + dц повторяет распределение dт, но максимальное значение dц на два порядка ниже, чем максимальное значение dт.
Для материала с относительно толстыми швами, низкопластичного по отношению к сдвигам (с высоким пределом т0) и высокопластичного по отношению к характеристикам кривизны (с низким пределом ц0), линии уровня диссипативной энергии dц приведены на рис. 5. Распределение этой величины сильно неоднородно и практически не зависит от толщины швов. На рис. 6 представлены поля моментных напряжений ц31 и ц32, которые также слабо зависят от 8. Масштаб неоднородности этих полей связан с внутренним параметром континуума — характерным масштабом микроструктуры материала.
Детальный анализ численных решений показывает, что моментные напряжения и связанные с ними кривизны слабо влияют на конечное макромасштабное деформированное состояние массива, которое характеризуется основными величинами — перемещениями и
кДж/м3 10.0 8.3 6.7 5.0 3.3 1.7 0.0
dт, кДж/м3 12 10
Рис. 4. Линии уровня энергии пластической диссипации для материалов с тонкими (а) и толстыми (6) швами
Ц32, кПа • м
■-5.5
-2.8 0.0 2.8 5.7 8.5
Рис. 6. Линии уровня моментных напряжений ц31 (а) и ц32 (6)
соответствующими им напряжениями. Распределение моментных напряжений имеет ячеистую структуру, отражающую неоднородность материала и изменение неоднородности в процессе нагружения. Поэтому, в отличие от обычных напряжений, их следует связывать с мезомасштабным уровнем деформации структурно-неоднородного материала. Хаотичное распределение энергии пластической диссипации за счет необратимого изменения тензора кривизны во всем объеме среды подтверждает гипотезу о том, что пластификация материала на мезоуровне обусловлена вращательными степенями свободы частиц.
Работа выполнена при финансовой поддержке Beijing High-caliber Talent from Overseas (BHT0 201612129-WD).
Литература
1. Панин B.E., Фомин B.M., Tumoв B.M. Физические принципы мезо-
механики поверхностных слоев и внутренних границ раздела в деформируемом твердом теле // Физ. мезомех. - 2003. - Т. б. -№ 2. - С. 5-14.
2. Панин B.E., Панин A.B., Mouœerno Д.Д. ««Шахматный» мезоэффект
интерфейса в гетерогенных средах в полях внешних воздействий // Физ. мезомех. - 2006. - Т. 9. - № б. - С. 5-15.
3. Панин B.E., Eгopyшкuн B.E., Панин A.B. Нелинейные волновые процессы в деформируемом твердом теле как многоуровневой иерархически организованной системе // УФН. - 2012. - Т. 182. -№ 12. - С. 1351-1357.
4. Панин B.E., Лuxaчeв B.A., Гpuняeв Ю.B. Структурные уровни деформации твердых тел. - Новосибирск: Наука, 1985. - 229 с.
5. Гу3ee M.A., Mampoe B.B. Деформирование и разрушение сильно сжатых горных пород вокруг выработок. - Владивосток: Даль-наука, 2007. - 231 с.
6. Rattez H., Stefanou I., Sulem J., Veveakis M., Poulet T. Numerical analysis of strain localization in rocks with thermo-hydro-mechanical couplings using Cosserat continuum // Rock Mech. Rock Eng. -2018. - V. 51. - No. 10. - P. 3295-3311. - doi 10.1007/ s00603-018-1529-7.
7. Tarasov B.G. Hitherto unknown shear rupture mechanism as a source of instability in intact hard rocks at highly confined compression // Tectonophysics. - 2014. - V. 621. - P. 69-84. - doi 10.1016/j.tecto. 2014.02.004.
8. Tarasov B.G., Guzev M.A., Sadovskii V.M., Cassidy M.J. Modelling the mechanical structure of extreme shear ruptures with friction approaching zero generated in brittle materials // Int. J. Fracture. - 2017. -V. 207. - No. 1. - P. 87-97. - doi 10.1007/s10704-017-0223-1.
9. Cosserat E., Cosserat F. Théorie des Corps Déformables. Chwolson's Traité Physique. - Paris: Librairie Scientifique A. Hermann et Fils, 1909. - P. 953-1173.
10. Gunther W. Zur Statik und Kinematik des Cosseratschen Kontinuums // Abhandlungen der Braunscheigschen Wissenschaftlichen Gesellschaft. - Braunschweig: F. Vieweg and Sohn, 1958. - V. 10. - P. 195213.
11. Sadovskii V.M., Sadovskaya O.V. Modeling of elastic waves in a blocky medium based on equations of the Cosserat continuum // Wave Motion. - 2015. - V. 52. - P. 138-150. - doi 10.1016/j.wavemoti. 2014.09.008.
12. Friedrichs K.O. Symmetric hyperbolic linear differential equations // Commun. Pure Appl. Math. - 1954. - V. 7. - No. 2. - P. 345-392. -doi 10.1002/cpa.3160070206.
13. Годунов C.K. Уравнения математической физики. - M.: Наука, 1979. - 391 с.
14. Godunov S.K., Romenskii E.I. Elements of Continuum Mechanics and Conservation Laws. - New York: Kluwer Academic / Plenum Publishers, 2003. - 260 p. - doi 10.1007/978-1-4757-5117-8.
15. Садовский B.M. Разрывные решения в задачах динамики упруго-пластических сред. - M.: Наука, 1997. - 208 с.
16. Sadovskaya O., Sadovskii V. Mathematical Modeling in Mechanics of Granular Materials. : Advanced Structured Materials, V. 21. - Heidelberg: Springer, 2012. - 390 p. - doi 10.1007/978-3-642-29053-4.
17. Свидетельство о гос. регистрации программы для ЭВМ № 2012614823 от 30.05.2012 (Роспатент). Программный комплекс для расчета двумерных динамических задач моментной теории упругости (2Dyn_Cosserat) / В.М. Садовский, О.В. Садовская, М.П. Варыгина. - М.: ФИПС, 2012. - RU ОБПБТ №3(80).-С. 209.
18. Тюменцев А.Н., Дитенберг И.А., Коротаев А.Д., Денисов К.И. Эволюция кривизны кристаллической решетки в металлических материалах на мезо- и наноструктурном уровнях пластической деформации // Физ. мезомех. - 2013. - Т. 16. - № 3. -С. 63-79.
Поступила в редакцию 08.02.2019 г., после доработки 18.02.2019 г., принята к публикации 18.02.2019 г.
Сведения об авторах
Садовский Владимир Михайлович, д.ф.-м.н., проф., дир. ИВМ СО РАН, проф. СФУ, sadov@icm.krasn.ru
Гузев Михаил Александрович, д.ф.-м.н., акад. РАН, проф., дир. ИПМ ДВО РАН, проф., Beijing University of Civil Engineering and Architecture, China, guzev@iam.dvo.ru
Садовская Оксана Викторовна, к.ф.-м.н., снс ИВМ СО РАН, o_sadov@icm.krasn.ru
Qi Chengzhi, Dr. Sci. (Phys. and Math.), Prof., Beijing University of Civil Engineering and Architecture, China, qichengzhi65@163.com
