CC BY
27
УДК 539.3
Автосолитонная концепция сейсмического процесса. Часть 2. Численные исследования генерации и распространения медленных деформационных автосолитонных возмущений
П.В. Макаров1,2, И.Ю. Смолин1,2, Ю.А. Хон2, М.О. Еремин1,2,
Р. А. Бакеев1,2, А.Ю. Перышкин1,2, В. А. Зимина1,2,
1 1 12 А. Чирков , А.А. Казакбаева , А.Ж. Ахметов '
1 Национальный исследовательский Томский государственный университет, Томск, 634050, Россия 2 Институт физики прочности и материаловедения СО РАН, Томск, 634055, Россия
В рамках автосолитонной концепции сейсмического процесса показано, что его полной математической моделью как процесса деформирования и разрушения нагружаемой геосреды являются динамические уравнения механики деформируемого твердого тела вместе со специфическими определяющими уравнениями, описывающими реологию геосреды. Эти уравнения описывают как обычную эволюцию напряженно-деформированного состояния, обусловленную распространением волн напряжений со скоростями звука, определяемыми спецификой определяющих уравнений, так и медленную динамику нагружаемой прочной среды. Численно изучены особенности генерации деформационных автосолитонов, структура фронтов и особенности распространения внутриразломных и межразлом-ных деформационных возмущений. Выполнено численное моделирование медленных деформационных возмущений в элементах реальных геосред.
Ключевые слова: автоволны, автосолитоны, неупругая деформация, землетрясения, моделирование, активная диссипативная среда, геосреда
DOI 10.24412/1683-805X-2021-1-18-36
Autosoliton view of the seismic process. Part 2. Possibility of generation and propagation of slow deformation autosoliton disturbances
in geomedia
P.V. Makarov1,2, I.Yu. Smolin1,2, Yu.A. Khon2, M.O. Eremin1,2, R.A. Bakeev1,2, A.Yu. Peryshkin1,2, V.A. Zimina1,2, A. Chirkov1, A.A. Kazakbaeva1, and A.Zh. Akhmetov1,2
1 National Research Tomsk State University, Tomsk, 634050, Russia 2 Institute of Strength Physics and Materials Science SB RAS, Tomsk, 634055, Russia
In the autosoliton view, the complete mathematical model of the seismic process taken as the deformation and fracture process of a loaded geomedium combines dynamic equations of solid mechanics and specific constitutive equations for geomedium rheology. These equations describe both the conventional stress-strain evolution due to the stress wave propagation with sound velocities, which are governed by special features of constitutive equations, and slow dynamics of the loaded strong medium. Numerical investigation is given to the generation of deformation autosolitons, front structure, and propagation of intra- and interfault deformation disturbances. Slow deformation disturbances in real geomedium elements are numerically modeled.
Keywords: autowaves, autosolitons, inelastic deformation, earthquakes, modeling, active dissipative medium, geomedium
© Макаров П.В., Смолин И.Ю., Хон Ю.А., Еремин М.О., Бакеев Р.А., Перышкин А.Ю., Зимина В.А., Чирков А., Казакбаева А.А., Ахметов А.Ж., 2021
1. Введение
В первой части статьи приведены обоснования для рассмотрения медленных деформационных возмущений в геосредах как автосолитонных решений соответствующих нелинейных математических моделей. Разные классы решений такого типа (бегущие, статические, пульсирующие) и их сопоставление с различными проявлениями деформационных явлений в земной коре позволяют сформулировать автосолитонную концепцию сейсмического процесса в целом. Мы показали, что медленные деформационные возмущения по своей физической природе являются бегущими деформационными автосолитонами. Также было показано, что разлом, точнее зона разлома, является стационарным автосолитоном (по современным представлениям разлом является геологическим телом — протяженной областью, обладающей существенной шириной; прочностные характеристики в этой области заметно меньше прилегающих консолидированных блочных структур). Существенной характеристикой разлома как стационарного автосолитона является то, что скорость деформации в теле разлома значительно выше, чем в прилегающих геоблоках. Также разлом как автоструктура является и пульсирующим автосолитоном [1, 2], когда скорость деформирования в отдельных его областях возрастает на многие порядки во время землетрясения. В такой активной среде, как нагруженная геосреда, возможно также развитие спонтанных процессов локальной потери устойчивости и образования новых автоструктур [1-3].
Для объяснения возможных свойств медленных деформационных возмущений в геологических средах предложены различные модели. Имеются модели медленной динамики, в которых уравнения движения постулируются. В частности, рассматриваются среды, описывающиеся уравнениями в частных производных параболического типа (диффузии) или уравнением синус-Гордона и его различными модификациями. Эти уравнения не связаны с деформационными процессами и эволюцией напряженно-деформированного состояния. По этой причине они не отражают физическую суть деформационных явлений в нагружаемой среде. Подобные модели являются имитационными, параметры в таких моделях подбираются по соответствующим данным наблюдений, свои для каждого конкретного случая.
Поэтому в этой части статьи будет представлена модель, основанная на полной системе уравне-
ний механики деформируемого твердого тела. Будет показано, что при определенных условиях, накладываемых на запись определяющих уравнений упругопластической среды, решения этих уравнений также содержат и медленные деформационные возмущения, которые распространяются на фоне сложившегося напряженно-деформированного состояния и упругих волн напряжений.
Поскольку надежных измерений медленных деформационных возмущений в геосредах крайне мало, в работе подробно численно изучен аналог таких медленных деформационных возмущений, наблюдаемых в нагружаемых металлических образцах, — формирование и распространение фронтов Людерса. Фронты Людерса по своей физической природе являются автоволнами, т.е. частным видом автосолитонов — простейшими бегущими деформационными автосолитонами. В численных расчетах будет подробно исследована структура фронта таких деформационных возмущений.
На примерах модельных геосред будут представлены особенности генерации, распространения и структуры фронтов медленных деформационных возмущений, распространяющихся как внутри разломов, так и внутри блоков от разломов, влияние силы тяжести на особенности деформационных фронтов в верхнем слое земной коры.
2. Зависимость скорости распространения фронта Людерса и длины площадки текучести от параметров математической модели
2.1. Введение
Закономерности распространения фронтов Лю-дерса хорошо изучены как экспериментально, так и теоретически. Однако не хватает работ, в которых подробно рассматривается феномен полос Людерса с точки зрения зависимости длины площадки текучести и скорости распространения полосы Людерса от параметров применяемых определяющих уравнений. В данной работе мы пытаемся восполнить пробел и провести комплексное численное исследование влияния параметров модели на отклик нагруженных образцов. Модели уп-ругопластического отклика среды, рассмотренные в [4-6], принципиально похожи. Их объединяет наличие зуба текучести и площадки текучести с двух-предельным критерием текучести, хотя и существуют некоторые специфические различия. Здесь мы объединяем идеи и уравнения этих моделей для проведения параметрического исследования.
Рис. 1. Типичная диаграмма пластического деформирования для низкоуглеродистой стали на начальных стадиях пластической деформации
Судя по макроскопической кривой течения а-в, зарождение полос локализованной пластической деформации происходит при достижении верхнего предела текучести с последующим формированием площадки текучести при относительно постоянном напряжении и при более низком пределе текучести [4, 7] (рис. 1).
В работе рассматриваются зависимости длины площадки текучести и скорости распространения фронта Людерса как деформационного автосоли-тонного возмущения от параметров математической модели, а именно от амплитуды сброса напряжения и коэффициента деформационного упрочнения.
По нашим представлениям, как процесс формирования медленных деформационных возмущений в нагруженной геосреде, так и процесс формирования фронтов Людерса в нагружаемых металлических образцах по своей физической природе являются процессами формирования деформационных автосолитонов в прочной нагружаемой среде. По этой причине фронты Людерса, условия и особенности их формирования и распространения являются корректной физической моделью медленных деформационных возмущений в геосредах. Численное моделирование и изучение зависимости их свойств от реологических особенностей упру-гопластической среды может существенно улучшить понимание формирования и распространения в геосредах подобных медленных деформационных автосолитонов. В этом контексте особое значение приобретает изучение процессов формирования полос локализованной деформации как стационарных автосолитонов — физических образов разломов.
Рис. 2. Схема образца, подвергнутого растягивающей нагрузке. Размеры рабочей части образца 50 х 10 х 2 мм
Изучение особенностей взаимодействия бегущих деформационных автосолитонных возмущений со стационарными деформационными автосо-литонами (разломами), слияние с ними и наращивание в них степени неупругой деформации (или поврежденности) может явиться ключом к пониманию процессов активизации разломов и перехода их эволюции в катастрофическую стадию.
2.2. Материалы и методы
В представленных численных экспериментах моделировался процесс деформирования малых образцов с размером рабочей части 50 х 10 х 2 мм3 и физико-механическими свойствами низкоуглеродистой стали (рис. 2, табл. 1). Для простоты моделирования мы рассматриваем только рабочую часть образца, в которой наблюдались фронты Людерса. Предполагается, что материал изотропен на макроуровне (рис. 2).
Моделирование нагружения проводилось в рамках трехмерного конечно-разностного метода. Он основан на явной схеме интегрирования по времени динамических уравнений механики твердого тела [8]. Схема устойчива при выполнении условия Куранта [8].
2.3. Математическая модель
Математическая модель, описывающая процессы неупругого деформирования и разрушения при наличии фронтов Людерса, представлена системой уравнений механики деформируемого твердого тела [8-12]. Это уравнения сохранения массы (1) и импульса (2). Система также содержит геометрические соотношения для тензоров полной скорости деформации (3) и завихренности (4). Среда считается баротропной. По этой причине система уравнений замкнута без привлечения закона сохранения энергии:
Таблица 1. Физико-механические свойства материала
р, г/см3 К, ГПа ц, ГПа
7.846 172 79.2
PV = PoVo, ■
ôxj
2s j = Vh j + v u, 2cc j = Vh j - Vu.
(1) (2)
(3)
(4)
ô t =é e +oP
f К, ôPq) = H/'
Y =
|cup, если epq = 0, |clow + если8Pq > 0,
up _low 5 F = 5
-Ас,
(8)
(9)
(10)
Здесь р0, р, У0, V — начальное и текущее значения плотности и объема соответственно; иг — компоненты вектора скорости; агу — компоненты тензора напряжений Коши; хг — декартовы координаты.
В качестве определяющих уравнений используются соотношения гипоупругого изотропного континуума. В рамках предположения об изотропной среде в уравнении состояния есть только две константы, а именно объемный К и сдвиговой ц модули. Тензор напряжений разбивается на шаровую (5) и девиаторную (6) части:
Р = -К 0,0 = в й, (5)
^у + ' - ' гк = У - 3 0 Ч' Р ), (6)
где Р — гидростатическое давление; 0 — объемная деформация; — тензор девиаторных напряжений; ер — тензор пластической деформации, точка над символами обозначает производную по времени. Коротационная производная Яуманна в уравнении (6) применяется для вычитания вращения ячейки (частицы или элемента среды) как целого.
Определяющие уравнения, описывающие пластические деформации, требуют обсуждения. Отметим, что мы предполагаем независимость отклика упругого континуума от пластического течения. В этом случае полная скорость деформации делится на две части (уравнение (7)), где индексы 1, е и р соответствуют полной, упругой и пластической деформации соответственно:
(7)
Как часть определяющего отклика используется двухпредельный критерий текучести, основанный на критерии текучести фон Мизеса (8), который широко применяется для металлов [6, 8]:
где f— функция текучести; spq — эквивалентная пластическая деформация; т — второй инвариант тензора девиаторных напряжений; Y — текущий предел текучести; cup — верхний предел текучести; clow — нижний предел текучести; k — коэффициент деформационного упрочнения.
Граничные и начальные условия:
- задаются скорости узлов на противоположных краях Z-плоскости: vz = v(t) и vz = —v(t) соответственно;
- все другие грани образца являются свободными поверхностями;
- ограничения скольжения узлов осуществляются нулевыми компонентами vx и vy вектора скорости на гранях, где заданы компоненты вектора скорости.
Исходное состояние образца соответствует нулевым значениям всех параметров напряженно-деформированного состояния.
При определении функции v(t) использовался прием медленного нагружения, чтобы минимизировать возможный динамический эффект. Примечательно, что время увеличения нагрузки соответствует более чем 20 пробегам упругой продольной волны через всю расчетную область, что дает практически квазистатическое условие эволюции напряжений и деформаций.
Приведенные ниже результаты представляют собой параметрическое исследование, в котором варьировались только два параметра модели, а именно: коэффициент деформационного упрочнения k в диапазоне от 25 до 300 МПа и разница между верхним и нижним пределом текучести Ас в диапазоне от нуля до 40 МПа при прочих равных параметрах (clow = 125 МПа).
2.4. Анализ результатов
В рассматриваемом случае изотропного материала две сопряженные полосы локализованной пластической деформации, наклоненные под углом «54°, зарождаются вблизи захватов и далее распространяются по образцу во всех выполненных расчетах (см. состояние 1 на рис. 3, б для примера). Стоит отметить, что в экспериментах не наблюдается симметричная картина крестообразных полос, которая представлена на рис. 3. Симметричная картина локализации пластической деформации получена в условиях полной симметричной постановки задачи, что свидетельствует о корректности численного решения и высокой точности расчета. Учет структурной неоднородности или
Рис. 3. Кривая напряжение-деформация для До = 40 МПа и коэффициента деформационного упрочнения к = 100 МПа (а); распространение полос Людерса в последовательные моменты 1-4 в модельном образце при одноосном растяжении (б)
отсутствие строгой симметрии в геометрии образца приводит к тому, что нарушается симметрия и в получаемом решении — вместо двух перекрещенных полос образуется одна от той угловой области, где раньше начались пластические деформации.
На рис. 3, а представлена диаграмма нагруже-ния образца. Зарождение полос Людерса вблизи обоих захватов происходит при достижении осевым растягивающим напряжением значения наблюдаемого верхнего предела текучести («214 МПа). Начало стадии площадки текучести характеризуется резким падением осевого напряжения в образце ниже значения нижнего предела текучести. Распространение полос Людерса происходит при относительно неизменном уровне осевого напряжения, соответствующего нижнему пределу текучести («190 МПа). При этом наблюдаются характерные осцилляции напряжения, связанные с неравномерным распространением полос локализованной пластической деформации и образованием волн разгрузки при релаксации напряжений в частицах материала, переходящих в пластическое состояние. По мере приближения к завершению стадии площадки текучести амплитуда осцилляций становится меньше, что свидетельствует об увеличении общего объема материала, перешедшего из упругого в пластическое состояние. Состояние 4
на рис. 3, б соответствует моменту нагружения, когда встречные фронты Людерса сталкиваются и аннигилируют. После этого на кривой течения начинается стадия параболического упрочнения.
Рассмотрим более подробно стадию распространения полос Людерса как волн переключения из упругого в пластическое состояние. Для этого проанализируем распределения интенсивности пластической деформации (рис. 4, а) и интенсивности скорости пластической деформации (рис. 4, б) вдоль линии АА', отмеченной на рис. 3, б (состояние 4). Распределения (состояния 1-4), представленные на рис. 4, соответствуют состояниям, представленным на рис. 3. За фронтом отчетливо виден след пластической деформации, значения которой слабо меняются вплоть до завершения стадии площадки текучести и в среднем составляют 0.012, за исключением всплесков на границах образца, обусловленных наличием концентраторов. Это означает, что вся деформация образца сосредоточена на фронте движущихся полос Людерса.
Согласно полученным распределениям, интенсивность скорости пластической деформации представляет собой возмущение, охватывающее приблизительно 10% от общей длины образца. При этом за пределами данного возмущения скорость пластической деформации отсутствует или близка к нулю. Продвижение фронта вдоль оси образца обеспечивается последовательными пробегами возмущения поперек образца от одной боковой грани к другой, что будет подробно обсуждено в разделе 3, приблизительно под углом в 54° в настоящем расчете. Пока фронт (его поперечная составляющая) не дошел до сечения АА', значения скорости неупругой деформации минимальны. Они максимальны при пересечении линии АА', затем опять уменьшаются, что и демонстрирует рис. 4, б (состояния 1 и 2).
Для выяснения характера зависимости длины площадки текучести, а также скорости распространения полос Людерса от параметров модели был выполнен численный эксперимент, в котором параметры модели Да и к варьировались с определенным шагом. На рис. 5 продемонстрированы кривые нагружения образцов при различных значениях сброса напряжения Да, но при фиксированном значении коэффициента деформационного упрочнения к. На основе данных моделирования можно отметить, что чем больше сброс напряжения Да, тем больше оказывается амплитуда зуба текучести и длина площадки текучести.
Рис. 4. Распределение интенсивности пластической деформации (а) и интенсивности скорости пластической деформации (б) по профилю АА для разных моментов времени 1-4
Рис. 5. Кривые напряжение-деформация для различных Да при постоянном значении коэффициента деформационного упрочнения к = 100 МПа, Да = 10 (1), 30 (2), 40 МПа (3)
Результаты параметрического исследования скорости распространения полосы Людерса как авто-солитонного возмущения и продолжительности площадки текучести в зависимости от величин Да и к показаны на рис. 6, а. Отношения скорости распространения фронта к скорости нагружения Рй/Рьаа для различных Да и к отражены на рис. 6, б. Следует отметить, что при малых значениях Да, не превышающих 5 МПа, для всех к соотношение между скоростью распространения полосы и скоростью нагружения остается достаточно высоким. Это свидетельствует о том, что волна переключения распространяется по всему образцу с высокой скоростью Уь = пУ\оа& (при п = 100-300). При дальнейшем увеличении Да наблюдается заметное изменение скорости распространения, приближающееся к наблюдаемому в экспериментальных исследованиях, например [13, 14]. Зависимости длины площадки текучести от Да являются нелинейно возрастающими функциями. Можно отметить, что чем ниже значение коэффициента де-
формационного упрочнения к, тем больше площадка текучести при фиксированном значении Да. Это свидетельствует о том, что продолжительность стадии площадки текучести зависит от скорости восстановления среднего уровня напряжений в нагруженном образце за фронтом пластической волны. Можно также отметить, что наблюдаемая в экспериментах разница между верхним и нижним пределами текучести в 2-2.5 раза меньше, чем истинный сброс напряжения, который реализуется в отдельных точках нагруженного образца.
3. Структура фронта деформационного автосолитонного возмущения
Для выявления сложной структуры деформационного автосолитона, распространяющегося вдоль образца, были численно исследованы особенности самоорганизованного распространения медленных деформационных возмущений неупругой природы в среде с упругопластическим переходом. Все численные расчеты проведены по модели, описанной в разделе 2.3.
Подробно проследить за изменением структуры фронта медленного деформационного автосолито-на можно по последовательности формирования очагов упругопластических переходов. Усредненная картина продвижения вдоль плоского образца деформационного возмущения — перехода из упругого состояния в пластическое — образует фронт автосолитонного возмущения. Локальные очаги такого перехода распространяются поперек образца в направлении максимальных касательных напряжений приблизительно под углом 45° к оси образца для случая плоского образца. Описанная структура фронта хорошо выявляется в скоростях неупругой деформации, как показано на рис. 7.
Да, МПа Да, МПа
Рис. 6. Кривые зависимости длины площадки текучести (а) и отношения скорости распространения полосы Людерса к скорости нагружения (б) от Да для разных значений к = 25 (1), 50 (2), 100 (3), 150 (4), 200 (5), 300 МПа (6)
Рис. 7. Структура фронта медленного деформационного автосолитона в модельном образце
Белые стрелки указывают на место в образце, где в данный момент находятся области локальных очагов, распространяющихся поперек образца. Черные стрелки указывают продольное направление, в котором распространяется медленный деформационный фронт.
На рис. 8 представлена схема последовательных этапов формирования и распространения фронта деформационного автосолитона. С течением времени локальное деформационное возмущение последовательно распространяется от одного края образца к другому под углом приблизительно 45°, формируя фронт медленного деформационного автосолитона и обеспечивая его продвижение вдоль образца. Локальные деформационные возмущения (уединенные импульсы скорости неупругой деформации) показаны жирными точками, у которых помещены символы, означающие последовательные моменты времени За фронтом ав-
Рис. 8. Схематическое представление структуры фронта медленного деформационного автосолитона с течением времени в образце
тосолитона остается область, охваченная пластической деформацией, она показана черным цветом.
Распространение поперек образца локальных возмущений скорости неупругой деформации происходит примерно на порядок быстрее распространения формируемого ими фронта медленного деформационного автосолитона вдоль образца.
Движущими силами деформационного автосо-литона являются максимальные сдвиговые напряжения, в направлении которых распространяется импульс скорости неупругих деформаций, пересекая образец от одной свободной границы до другой. Это можно наблюдать на рис. 9. Здесь представлены совмещенные картины распределения главных сдвиговых напряжений и интенсивности скоростей неупругой деформации. Наклонные стрелки указывают на изолинии максимальных касательных напряжений, куда в последующий момент времени продвинется очаг упругопластического перехода. При этом наблюдается не только движение очага возмущения скорости пластической деформации, но и движение наклонного фронта медленного деформационного автосолитона вдоль оси образца (оси приложения нагрузки), как показано горизонтальной стрелкой у верхней грани образца на рис. 9.
Таким образом, распространение деформационного автосолитонного возмущения вдоль оси об-
Рис. 9. Иллюстрация структуры фронта медленного деформационного автосолитона на этапе его распространения (цветной в онлайн-версии)
Рис. 10. Динамика фронта неупругой деформации в узком пологом разломе при сжатии геосреды для последовательных моментов времени: 1, 2 — прочные блоки, 3 — узкий пологий разлом. Черным цветом отмечена область разлома, охваченная пластической деформацией
разца в направлении нагружения (сжатие либо растяжение) обеспечивается последовательным прохождением возмущения скорости неупругой деформации поперек образца в направлении максимальных касательных напряжений (под углом 45° в данном случае плоского образца). За фронтом остается след в виде накопленной неупругой деформации, величина которого определяется амплитудой нагружения.
4. Моделирование внутриразломных автосолитонных возмущений
Задача распространения деформационного возмущения вдоль разломной зоны решалась численными методами механики деформируемого твердого тела (уравнения (1)-(8)) совместно с дискретным методом клеточных автоматов, который подробно описан в работах [9, 15]. По этой методике каждая расчетная ячейка (частица среды) переходит из упругого состояния в пластическое, порождая деформационный автосолитон, только тогда, когда напряжения в ней достигнут верхнего предела и хотя бы в одной из соседних ячеек накопленная неупругая деформация также достигнет заданного критического значения.
Для изучения процесса распространения медленного автосолитонного возмущения по разлому в стесненных условиях были проведены две группы численных экспериментов. В первой из них нагружаемая модель представлена двумя блоками, разделенными наклонным разломом малой толщины (рис. 10). На рис. 10 приведена динамика фронта неупругой деформации в пологом разломе при сжатии геосреды. Нагрузка с постоянной скоростью прилагалась на верхней грани, смещения боковых граней по горизонтали запрещены. Видно, что фронт неупругой деформации распространяется от точки с наибольшей концентрацией напряже-
ния (левый нижний угол). Наблюдается отчетливо выраженный предвестник, бегущий впереди возмущения по нижней границе разлома.
Во второй серии численных экспериментов нагружаемая среда представлена двумя граничными блоками с находящимся между ними достаточно широким разломом, ориентированным под малым углом к оси приложения нагрузки (рис. 11).
Моделируемая среда представляет собой протяженную область разлома, обладающего пониженными прочностными и упругими характеристиками и окружающей его более прочной средой. Прочность геосреды в разломной зоне задана в 3 раза ниже окружающей среды.
Боковые границы могут перемещаться в вертикальном направлении, а нижняя граница — в горизонтальном направлении. На верхней границе задано сжимающее усилие модельной среды с постоянной скоростью до момента времени в 700 мкс, порождая импульс сжатия, необходимый для генерации деформационного возмущения, далее скорость ее движения уменьшается в 10 раз для предотвращения перегрузки среды. Это уменьшенное воздействие необходимо для поддержания им-
Рис. 11. Схема расчетной области — широкий разлом с двумя прилегающими блоками
Рис. 12. Деформационный фронт в широком разломе (показана только область разлома). Картины распределения накопленной деформации (а) и скорости неупругой деформации (б) в последовательные моменты времени (значение времени указано в условных единицах): 11 = 678, ^ = 956, Ц = 1140, t4 = 1260, t5 = 1319. Стрелками указаны направления распространения отдельных фронтов
пульса деформационного возмущения. Эти условия имитируют возможный реальный процесс генерации медленных автосолитонных возмущений в геосредах. На фоне медленного нагружения блочной геосреды возникает краткий импульс. В настоящих расчетах он в 10 раз выше последующего фонового сжатия.
В широкой разломной зоне распространение фронта деформационного возмущения принципиально отличается от динамики фронта в узком разломе. Его распространение вдоль разлома существенно нестационарно, так на первых этапах в областях интерфейсов (граница между разломом и прилегающими блоками на рис. 12) формируются узкие зоны с повышенными значениями скорости неупругой деформации (рис. 12, б). Эти области ориентированы в направлении максимальных касательных напряжений. Они разрастаются поперек разлома в виде треугольных областей с предвестниками, движущимися вдоль разлома, в областях его сопряжения с прочными блоками (вдоль интерфейсов). Характерно, что в разломной зоне остаются области, практически неохваченные упру-гопластическим переходом (светлые участки на рис. 12, а). Такие слабо деформированные области обусловлены релаксацией напряжений за разрастающимися областями неупругой деформации. Таким образом, пока треугольные области неупруго деформированной среды, разрастаясь, пересекают область разлома, продвижение фронта деформаци-
онного возмущения вдоль разлома незначительно. Многократное повторение описанного процесса разрастания треугольных областей, охваченных упругопластическими переходами, обеспечивает распространение возмущения вдоль такого широкого разлома. Зарождение возмущения всегда происходит на интерфейсах. В приведенном примере распространение возмущения обеспечивается подкачкой энергии в среду — фоновым сжатием на верхней границе, а также стесненной деформацией (запрет горизонтальных перемещений на боковых границах, прилегающих к разлому блоков), что и обеспечивает первоначальное зарождение возмущения на интерфейсах.
5. Численное моделирование процессов генераций и распространения деформационных межразломных автосолитонов
Для описания медленных фронтов деформации неупругой природы между разломами также использована модель с критерием текучести Друке-ра-Прагера, скомбинированная с методом клеточных автоматов [9, 15]. В качестве областей, вблизи которых могут зародиться фронты пластической деформации, заданы зоны разломов — узкие (тонкие) вытянутые области, имеющие ослабленные механические характеристики. Исследовалось поведение участка среды с разломами при одноосном сжатии вдоль вертикальной оси в условиях
Рис. 13. Распределение интенсивности неупругих деформаций при распространении деформационных возмущений в среде с разломом, ориентированным под углом 90° к оси сжатия, при реализации метода клеточных автоматов с окрестностями фон Неймана (а-в) и Мура (г-е)
плоской деформации при разной ориентации разломов: перпендикулярно (рис. 13) и параллельно оси сжатия (рис. 14).
Возможны разные реализации алгоритма клеточных автоматов для передачи медленного возмущения. Для плоских задач были использованы два варианта задания окрестности в двумерных моделях клеточных автоматов: фон Неймана и Мура. В окрестности фон Неймана соседями данной клетки являются клетки, имеющие с заданной общую сторону (общее ребро). Таким образом, каждая клетка имеет ровно четырех соседей. В окрестности Мура две клетки являются соседними, если они имеют либо общее ребро, либо общую вершину, т.е. в локальной окрестности каждой клетки (помимо ее самой) находится восемь клеток.
В связи с вышесказанным было проанализировано, как влияют на распространение деформационных автосолитонов между разломами их ориентация по отношению к оси нагружения и выбор варианта реализации алгоритмов клеточных автоматов.
На рис. 13 представлены картины распределения интенсивностей неупругих деформаций, отражающие распространение фронтов деформационных возмущений в нелинейных средах с разными вариантами реализации алгоритмов клеточных автоматов в случае, когда разлом ориентирован перпендикулярно направлению сжатия. Разломы на всех рисунках показаны в виде линий АВ.
Из рис. 13 видно, что фронты деформационных возмущений первоначально зарождаются на краях разломов, в областях концентраторов напряжений, имеют симметричную ромбовидную форму, которая несколько отличается в зависимости от выбранного алгоритма клеточных автоматов. Симметрия фронтов нарушается при их взаимодействии (рис. 13, б, в и д, е). Особо заметно отличие в случае перекрытия фронтов. В области перекрытия фронтов образуются множественные полосы локализации деформации. Скорость распространения медленных деформационных возмущений в среде с клеточными автоматами при заданных параметрах модели равна 25 и 125 м/с для окрестностей с четырьмя и восемью ячейками соответствен-
ООО 0.03 0.06 0.09 ООО 0.03 0.06 0.09 0 00 0.09 0.18 0.29
I ' ер, % I ' ер, % I ' ер, %
Рис. 14. Распределение интенсивности неупругих деформаций при распространении деформационных возмущений в среде с разломом, ориентированном под углом 0° к оси сжатия, при реализации метода клеточных автоматов с окрестностями фон Неймана (а-в) и Мура (г-е)
но. Как и следовало ожидать, для окрестности с восемью ячейками скорость распространения деформационных возмущений существенно выше.
При ориентации разлома вдоль оси сжатия (рис. 14) фронты сохраняют практически симметричную форму. При их взаимодействии, в области перекрытия также образуются полосы локализованной деформации, но они выражены слабее. Наибольшая деформация сосредоточена в двух полосах, пересекающихся в центре разлома (рис. 14, в, е).
Хронограммы распространения медленного деформационного автосолитона от одного наклонного разлома представлены на рис. 15. Так же, как и для ортогонального расположения разломов, в вершинах разлома генерируются медленные деформационные автосолитоны. Форма фронтов не ромбовидная, а треугольная, и они распространяются в противоположных направлениях к границам расчетной области. Еще больше отличие фронтов от наклонных разломов проявляется на поздних стадиях. Фронт автосолитона в вертикальном направлении (параллельно оси нагружения) распространяется быстрее, чем в горизонтальном на-
правлении. Более быстрая вертикальная вершина фронта является зародышем новых фронтов. Ввиду сдвигового характера пластических деформаций, фронт является не прямым, а зигзагообразным. Он «лавирует» в поле действия касательных напряжений. В зависимости от реализованного варианта окрестности клеточных автоматов, количество зигзагов меняется. Расчеты с окрестностью Мура приводят к образованию одного крупного треугольного фронта от каждого края разлома, они поглощают другие мелкие фронты. В случае окрестности фон Неймана количество треугольных фронтов больше и они разные по размеру. Иногда более крупные фронты поглощают мелкие.
В случае двух разломов в среде, ориентированных вдоль оси нагружения, картина взаимодействия фронтов представлена на рис. 16 и 17. Наличие другого разлома приводит к нарушению симметрии фронта еще до того, как фронты встретились. Это обусловлено (как и в предыдущих случаях) взаимодействиями звуковых волн напряжений, скорости которых на несколько порядков выше скоростей медленных деформационных фронтов. Эти взаимодействия формируют слож-
Рис. 15. Хронограммы распространения деформационных автосолитонов от разлома, наклоненного под углом 45° к оси нагружения, при реализации метода клеточных автоматов с окрестностями фон Неймана (а, б) и Мура (в, г). Показан след за фронтом в виде неупругой деформации
ную картину напряженно-деформированного состояния, влияющую на скорости и интенсивности деформаций в деформационных автосолитонных возмущениях.
На рис. 18 представлена хронограмма распространения медленных деформационных автосолито-нов, сформировавшихся на границах протяженных разломов АВ (рис. 18, а). Видно, что фронты де-
формационных возмущений приближаются друг к другу примерно с одинаковыми скоростями, их форма меняется с течением времени и при встрече они взаимодействуют, уменьшая область, не охваченную неупругой деформацией.
На рис. 13-18 показан след за фронтами деформационных автосолитонов в виде накопленной неупругой деформации.
Рис. 16. Хронограмма распространения деформационных автосолитонов, когда разломы расположены под углом 0° (для окрестности фон Неймана)
Рис. 17. Хронограмма распространения деформационных автосолитонов, когда разломы наклонены под углом 45° (окрестность фон Неймана). Приведены картины распределения неупругой деформации за фронтом
Рис. 18. Хронограмма распространения деформационных автосолитонов, когда разломы расположены вблизи границ среды
6. Деформационные фронты в геосредах
Проведено численное моделирование распространения в верхнем слое земной коры фронта деформационной волны, возникшего вследствие коллизии. Использована модель среды Друкера-Пра-гера с неассоциированным законом течения [1619]. Объект моделирования — верхний слой осадочного чехла толщиной 1.6 км. Расчет проведен на детальной сетке с шагом 10 м для точного определения формы деформационного фронта. Физико-механические свойства модельной среды приняты близкими к свойствам ослабленных горных пород (табл. 2).
В данном случае волна переключения упруго-пластического состояния описана с помощью двух-предельной кривой течения [20]. На рис. 19 представлен общий вид кривой. По горизонтали отложены абсолютные значения накопленной пластической деформации, по вертикали — отношение текущего значения когезии к начальному. Характерной особенностью этой зависимости когезии от
накопленной неупругой деформации является то, что после достижения предела упругости в локальной области происходит резкий сброс напряжений. Распространение пластической деформации на некоторое время приостанавливается, пока напряжения в этой области с ростом пластической деформации не вернутся к начальному значению когезии на ветви упрочнения. Подобный механизм распространения деформационного фронта реализован в модели с клеточными автоматами с той разницей, что в автоматах используется деформа-
Таблица 2. Физико-механические свойства песчаника
Плотность р, г/см3 2.2
Модуль объемного сжатия К, ГПа 12.8
Модуль сдвига О, ГПа 5.34
Начальная когезия У0, МПа 4
Коэффициент внутреннего трения а 0.1
Коэффициент дилатансии Р 0.01
т0
0.6-1—.—.—.—.—.—.—.—.—.— 0.000 0.004 0.008
Пластическая деформация, ед.
Рис. 19. Зависимость относительной когезии от накопленной пластической деформации
ционный критерий — переход в пластическое состояние происходит в том случае, когда у любого соседа пластическая деформация достигает критического значения.
На первом этапе моделирования проведен расчет установления гравитационных напряжений под действием силы тяжести в слое геосреды размером 1.6 х 10 км. На всех гранях, кроме верхней свободной поверхности, были заданы нулевые нормальные скорости. Параметры материала подобраны таким образом, что под действием гравитационных напряжений в моделируемом слое толщиной 1.6 км не происходит переход в пластическое состояние. Для этого требуется дополнительная нагрузка — тектоническое сжатие или растяжение.
На втором этапе моделирования на боковой границе задается нормальная скорость (достаточно низкая, чтобы избежать динамических эффектов). Это моделирует коллизионное взаимодействие тектонических плит. Максимальное значение скорости 1 м/с. Ограничения возможностей машинного счета требуют приема «сжатия времени» путем увеличения на несколько порядков скоростей на-гружения. Понятно, что на эти порядки увеличатся и скорости сгенерированных фронтов. При заданных скоростях сжатия скорость движения деформационного фронта достигает значений 700 м/с.
На рис. 20 показана динамика фронта. Хорошо видно, что фронт деформационной волны имеет форму «лезвия ножа». Объяснить это можно тем, что предел текучести по модели Друкера-Прагера растет с глубиной, а значит, более глубокие слои среды переходят в пластическое состояние с некоторой задержкой, связанной с тем, что уровень напряжений еще не достиг критического уровня. При этом важно отметить, что есть некоторая предельная глубина, ниже которой пластическое состояние в отсутствии тектонических напряжений не распространяется. В дальнейшем пластическая
Рис. 20. Распространение фронта деформационной волны по поверхности слоя геосреды (показано распределение скорости неупругой деформации)
деформация за счет сжатия распространится на всю глубину слоя геосреды. Но на начальном этапе коллизионного взаимодействия, фронт деформационной волны бежит только в приповерхностных слоях.
Для анализа механизмов формирования медленных движений было проведено численное моделирование процесса генерации медленных деформационных волн на интерфейсах блочной геологической среды в условиях быстрого смещения бортов разломов. На основе метода клеточных автоматов получено численное решение задачи распространения фронта пластической деформации в геосреде, находящейся под действием гравитационных напряжений, и в условиях действия на ее границе тектонических растягивающих либо сжимающих напряжений.
В качестве исследуемого объекта взята подобласть глубиной 10 км и протяженностью 50 км. Физико-механические свойства геосреды были выбраны на основе результатов глубинного сейсмического зондирования по одному из участков геологического профиля Батолит-1982 (табл. 3).
На первом этапе решена задача установления равновесного напряженного состояния под дейст-
Таблица 3. Физико-механические свойства геосреды
Плотность р, г/см3 2.7
Модуль сдвига О, ГПа 43
Модуль объемного сжатия К, ГПа 60
Предел прочности, МПа 10
Коэффициент внутреннего трения а 0.34
Коэффициент дилатансии р 0.087
вием силы тяжести. Затем область была подвергнута растяжению или сжатию через задание кинематических граничных условий на левой границе. Правая граница оставалась неподвижной.
На рис. 21 показано распространение фронта пластической деформации в слое геосреды при тектоническом растяжении. В случае дополнительного тектонического растяжения фронт зарождается у нижнего левого края горного массива (рис. 21, а), достигает дневной поверхности, затем распространяется вдоль горизонтальной оси и по всей глубине слоя. Такой характер формирования фронта обусловлен тем, что интенсивность касательных напряжений на глубине максимальная, а тектоническое растяжение вызывает разгрузку — снижение горизонтального сжимающего напряжения. Если до тектонической нагрузки эффективная прочность превышала гравитационные напряжения, то в результате разгрузки прочность снижается из-за падения давления и становится ниже интенсивности касательных напряжений. С увеличением смещения левой границы фронт пластической деформации достигает земной поверхности (рис. 21, б) и затем, сохраняя наклон, поступательно распространяется по слою вправо (рис. 21, в, г). Скорость движения фронта деформации сохраняется постоянной в течение всего процесса нагру-жения и составляет в данном случае 528 м/с.
Анализируя картины распространения фронта скорости пластической деформации (рис. 22), необходимо отметить свойственную модели Друке-ра-Прагера блочную (полосовую) картину развития пластической деформации в среде, как и на рис. 21. Формируются две системы сопряженных наклонных полос локализации пластической деформации, отстоящих друг от друга на некотором расстоянии. На картине распределения пластической деформации они ярко выражены и по мере распространения фронта пластической деформации вглубь среды остаются неподвижны. В них только увеличивается степень деформации. Фронт скорости пластической деформации движется также с постоянной скоростью, но за ним остаются отдельные статические зоны, соответствующие растущим полосам локализации. Кроме этого линия фронта может менять наклон в связи с изменением напряженного состояния (рис. 22, г).
В отличие от тектонического растяжения при тектоническом сжатии пластическая деформация зарождается у верхнего левого края горного массива (рис. 23). Это обусловлено тем, что рост горизонтального сжимающего напряжения приводит к росту давления и, как следствие, к росту эффективной прочности. При одновременном росте интенсивности касательных напряжений критическое состояние будет раньше достигнуто в верхних
Рис. 21. Распределение пластической деформации за фронтом пластической деформации при тектоническом растяжении в различные моменты времени: 70 (а), 120 (б), 295 (в), 560 с (г)
Рис. 22. Распространение фронта скорости пластической деформации при растяжении в различные момен ты времени: 70 (а), 120 (б), 295 (в), 560 с (г)
Рис. 23. Распространение фронта пластической деформации при тектоническом сжатии в различные моменты времени: 165 (а), 195 (б), 215 с (в)
Рис. 24. Распространение фронта скорости пластической деформации при сжатии в различные моменты времени: 165 (а), 195 (б), 215 с (в)
слоях геосреды. Фронт пластической деформации распространяется также вдоль слоя, но не на всю глубину, образуя ступенчатую область локализаций пластических деформаций, простирающуюся от земной поверхности до глубины 5-6 км. Скорость движения фронта 688 м/с. После достижения правой неподвижной границы структурной модели фронт пластической деформации начинает распространяться вглубь геосреды (рис. 23, в) со скоростью на порядок ниже — 67 м/с.
На рис. 24 представлено распределение скорости пластической деформации при сжатии слоя. Также можно выделить два этапа распространения фронта. На первом этапе фронт скорости пластической деформации распространяется в верхней части слоя геосреды в продольном направлении до глубины 5-6 км. На втором этапе при достижении фронтом правой неподвижной части горного массива фронт меняет направление и начинает распространяться вертикально вглубь горного массива (рис. 24, в).
7. Заключение
Наиболее общей математической моделью деформационных автосолитонных возмущений являются уравнения механики деформируемого твердого тела. Показано, что решения этих уравнений на больших временах могут приводить к генерации деформационных возмущений, которые с малыми скоростями распространяются на фоне сложившегося в прочной среде напряженно-деформированного состояния. Эти возмущения, как пра-
вило, низкоамплитудны. Как показало численное моделирование, скорости деформационных ав-тосолитонных возмущений пропорциональны скоростям импульсных воздействий на границах, сгенерировавших эти деформационные автосолитон-ные возмущения. Это свойство расчетных деформационных автосолитонных возмущений отвечает свойствам классических автосолитонов [1].
На примере моделирования деформационных фронтов Людерса, которые нами рассматриваются как один из видов типичных бегущих автосоли-тонных возмущений — автоволны, численно изучены зависимости параметров возмущений (скорости, амплитуды) от реологии прочной среды (верхнего и нижнего предела текучести, коэффициента деформационного упрочнения). Изучена структура фронтов, зависимости скорости фронта от скорости нагружения и реологии среды.
Подробно изучены особенности генерации, распространения и структуры фронтов как внутрираз-ломных, так и межразломных деформационных автосолитонных возмущений. В случае узких раз-ломных зон (для которых не выявляется структура фронта) внутриразломное автосолитонное возмущение распространяется по разлому со скоростью, пропорциональной скорости импульсного воздействия. Как показали расчеты, при прекращении этого воздействия скорость возмущения быстро падает, возмущение затухает.
В широких разломах формируются медленные автосолитонные возмущения сложной структуры. Изучена структура фронта такого возмущения при продольном сжатии среды, состоящей из двух бло-
ков с разломом между ними. Возмущение формируется в области концентраторов напряжений на интерфейсах и, разрастаясь, пересекает разлом. Эти процессы формирования медленных возмущений происходят последовательно на двух противоположных интерфейсах разлома, постепенно продвигаясь по разлому в продольном направлении.
Межразломные автосолитонные возмущения также первоначально формируются в областях концентраторов напряжений на интерфейсах разлом-блок. Численно изучены особенности формирования и распространения автосолитонных возмущений в блочной среде при различной ориентации разломов по отношению к направлению нагруже-ния — разлом перпендикулярен, параллелен, находится под углом к оси нагружения. Во всех этих случаях формируется сложная картина взаимодействующих деформационных автосолитонов. Неупругая деформация за фронтами существенно неоднородна, локализована. Области локализации неупругой деформации и повреждений можно рас-смативать как зародыши разломов разных масштабов.
Таким образом, изучение параметров и свойств медленных деформационных возмущений выявило их автоволновую природу. Полученные решения показали, что в нагружаемой прочной среде формируются три типа деформационных автосо-литонов: статические, пульсирующие, бегущие.
При определенных условиях бегущие автосоли-тоны взаимодействуют как не вполне упругие частицы, отталкиваются либо сливаются со статическими, формируя в прочной упругопластической среде диссипативные структуры — полосы локализованной деформации. Этот процесс приближает отдельные области разлома к критическому состоянию.
Полученные результаты позволили сформулировать автосолитонную концепцию медленных деформационных возмущений в Земле и процессов формирования в геосредах диссипативных структур, к которым следует отнести разломы. Разлом как статический автосолитон превращается в пульсирующий при активизации в нем или в его части деформационного процесса. Такие активизации разных масштабов — разных интенсивностей скоростей деформирования — можно трактовать как сейсмический процесс.
Работа выполнена при поддержке Российского научного фонда (проект № 19-17-00122).
Литература
1. Кернер Б.С., Осипов В.В. Автосолитоны // Успехи физических наук. - 1989. - Т. 157. - № 2. - С. 201261.
2. Васильев В.А., Романовский Ю.М., Яхно В.Г. Автоволновые процессы в распределенных кинетических системах // Успехи физических наук. - 1979. -Т. 128. - № 8. - С. 625-626.
3. Masterov A.V., Tolkov V.N., Yakhno V.G. Nonlinear Waves: Dynamics and Evolution. - New York: Sprin-ger-Velag, 1988.
4. Schwab R., Ruff V. On the nature of the yield point phenomenon // Acta Mater. - 2013. - V. 61. - No. 5. -P. 1798-1808.
5. Hallai J., Kyriakides S. On the effect of Luders bands on the bending of steel tubes. Part II: Analysis // Int. J. Solids Struct. - 2011. - V. 48. - No. 24. - P. 3285-3298.
6. Romanova V., Balokhonov R., Schmauder S. Three-dimensional analysis of mesoscale deformation phenomena in welded low-carbon steel // Mater. Sci. Eng. A. - 2011. - V. 528. - No. 15. - P. 5271-5277.
7. Johnson D., Edwards M., Chard-Tuckey P. Microstructural effects on the magnitude of Luders strains in a low alloy steel // Mater. Sci. Eng. A. - 2015. - V. 625. -P. 36-45.
8. Wilkins M. Computer Simulation of Dynamic Phenomena. - New York: Springer-Verlag, 1999.
9. Makarov P.V., Peryshkin A. Yu. Slow motions as inelastic strain autowaves in ductile and brittle media // Phys. Mesomech. - 2017. - V. 20. - No. 2. - P. 209-221. -doi 10.1134/S1029959917020114
10. Balokhonov R., Zinoviev A., Romanova V., Zinovie-va O. The computational micromechanics of materials with porous ceramic coatings // Meccanica. - 2016. -V. 51. - No. 2. - P. 415-428.
11. Smolin I., Makarov P., Eremin M., Matyko K. Numerical simulation of mesomechanical behavior of porous brittle materials // Proc. Struct. Integr. - 2016. - V. 2. -P. 3353-3360.
12. Eremin M. Numerical simulation of failure of sandstone specimens utilizing the finite-difference continuous damage mechanics approach // Proc. Struct. Integr. - 2019. - V. 18. - P. 135-141.
13. Danilov V., Gorbatenko V., Zuev L., Orlova D. Kinetics and morphology of Luders deformation in specimens with homogeneous structure and with a weld joint // Mater. Sci. Eng. A. - 2018. - V. 714. - P. 160-166.
14. Barannikova S., Ponomareva A.V., Zuev L.B., Veki-lov Yu.Kh., Abrikosov I.A. Significant correlation between macroscopic and microscopic parameters for the description of localized plastic flow auto-waves in deforming alloys // Solid State Commun. - 2012. -V. 152. - No. 9. - P. 784-787.
15. Makarov P.V., Khon Yu.A., Peryshkin A.Yu. Slow deformation fronts: model and features of distribution // Geodyn. Tectonophys. - 2018. -V. 9. - No. 3. -P. 755-769. - doi 10.5800/GT-2018-9-3-0370
16. Николаевский В.Н. Определяющие уравнения пластического деформирования сыпучей среды // Прикладная математика и механика. - 1971. - Т. 35. -№ 6. - С. 1017-1029.
17. Друкер Д., Прагер В. Механика грунтов и пластический анализ или предельное проектирование // Механика. Новое в зарубежной науке. Вып. 2. Определяющие законы механики грунтов. - М.: Мир, 1975. - С. 166-177.
18. Stefanov Yu.P., Bakeev R.A. Formation of flower structures in a geological layer at a strike slip displacement in the basement // Izv. Phys. Solid Earth. - 2015. -V. 51. - No. 4. - P. 535-547.
19. Stefanov Yu.P., Tataurova A.A. Effect of friction and strength properties of the medium on shear band formation in thrust structures // Phys. Mesomech. - 2019. -V. 22. - No. 6. - P. 463-472. - doi 10.1134/S10299599 19060031
20. Balokhonov R.R., Romanova V.A., Martynov S.A., Schwab E.A. Simulation of deformation and fracture of coated material with account for propagation of a Lu-ders-Chernov band in the steel substrate // Phys. Mesomech. - 2013. - V. 16. - No. 2. - P. 133-140.
Поступила в редакцию 14.12.2020 г., после доработки 12.02.2021 г., принята к публикации 16.02.2021 г.
Сведения об авторах
Макаров Павел Васильевич, д.ф.-м.н., проф., проф. ТГУ, гнс ИФПМ СО РАН, pvm@ispms.ru Смолин Игорь Юрьевич, д.ф.-м.н., доц., проф. ТГУ, зав. лаб. ИФПМ СО РАН, smolin@ispms.ru Хон Юрий Андреевич, д.ф.-м.н., проф., зав. лаб. ИФПМ СО РАН, khon@ispms.ru Еремин Михаил Олегович, к.ф.-м.н., снс ТГУ, нс ИФПМ СО РАН, eremin@ispms.ru Бакеев Рустам Альфредович, к.ф.-м.н., снс ТГУ, нс ИФПМ СО РАН, bakeev@ispms.ru Перышкин Алексей Юрьевич, мнс ТГУ, инж. ИФПМ СО РАН, alexb700@yandex.ru Зимина Валентина Алексеевна, мнс ТГУ, мнс ИФПМ СО РАН, miva@ispms.ru Чирков Артем, мнс ТГУ, chirkovartyem@gmail.com Казакбаева Айгерим Азаматовна, мнс ТГУ, aigerim_@bk.ru
Ахметов Аян Жанатович, мнс ТГУ, мнс ИФПМ СО РАН, ayan.akhmetov93@gmail.com
