CC BY
72
УДК 621.891
Возникновение поворотных мод деформации в приграничных зонах зеренной структуры в нагруженном поликристалле
Д.Д. Моисеенко, Ю.И. Почивалов, П.В. Максимов, В.Е. Панин
Институт физики прочности и материаловедения СО РАН, Томск, 634021, Россия
Проведены экспериментальные исследования формирования дефектной структуры в тонком поверхностном слое стали 20Х18Н10Т в процессе изнашивания в паре трения. Показано, что пластическая деформация в стесненных условиях пары трения является неоднородной и носит вихревой ротационный характер. Для теоретического описания процессов формирования вихревых структур в многоуровневых поликристаллических структурах развит метод возбудимых клеточных автоматов. Численные эксперименты показали формирование ротационных мод деформации на различных структурно-масштабных уровнях — модуляцию значений локального момента силы вдоль межзеренных границ и фрагментацию исходных структурных элементов в исходных конгломератах самосогласованно деформирующихся зерен в поле моментных напряжений.
Ключевые слова: физическая мезомеханика, трибоконтакт, компьютерное моделирование, клеточные автоматы
Rotational deformation modes in near-boundary regions of grain structure
in a loaded polycrystal
D.D. Moiseenko, Yu.I. Pochivalov, P.V. Maksimov, and V.E. Panin
Institute of Strength Physics and Materials Science, SB RAS, Tomsk, 634021, Russia
Experiments were performed to study the defect structure formation in a thin surface layer of 20Cr18Ni10Ti steel in frictional wear. It is shown that plastic deformation constrained by a friction pair is inhomogeneous and displays a rotational vortex character. A method of excitable cellular automata is developed for theoretical description of the formation of vortex structures in multiscale polycrystalline structures. Numerical experiments demonstrate the formation of rotational deformation modes on different structural scales — modulation of local couple forces along the grain boundary and fragmentation of initial structural elements in initial conglomerates of self-consistently deformed grains in the field of couple stresses.
Keywords: physical mesomechanics, tribocontact, computer simulation, cellular automata
1. Введение
Твердое тело как многоуровневая иерархическая система в процессе экстремальных внешних воздействий (механических, тепловых, радиационных и пр.) на каждом масштабном уровне претерпевает структурно-фазовые превращения, происходящие в результате действия потоков различной природы. Эти потоки можно разделить на два класса: потоки субстанции (частиц материи) и потоки энергии. К первому классу относятся потоки, осуществляющие массоперенос (диффузия, конвективно-адвективный перенос и пр.). Второй класс потоков является более широким и разнородным: это
могут быть перенос механической энергии в акустической волне, теплоперенос, электрический ток и др.
В работе [1] показано, что ранее теоретически предсказанные [2-4] и экспериментально обнаруженные [59] локализованные волны пластической деформации являются результатом сложной самоорганизации потоков энергии и локальных структурных трансформаций в иерархии структурно-масштабных уровней.
Следует обратить особое внимание на то, что, согласно современному квантовомеханическому представлению о материи, каждая элементарная частица является носителем кванта энергии. Поэтому вышеизло-
© Моисеенко Д.Д., Почивалов Ю.И., Максимов П.В., Панин В.Е., 2013
женное разделение потоков по типам является весьма условным, а выбор типов, определяющих тот или иной процесс потоков, зависит от специфики поставленной задачи и масштабного уровня процесса.
Проблема построения адекватной математической модели процесса переноса, которая учитывает иерархию масштабных уровней в структурно-фазовых превращениях, является актуальной задачей теории твердого тела. Подобный подход развивается в физической мезомеха-нике [10], в которой деформация нагруженного твердого тела рассматривается как многоуровневый процесс, связанный с потерей сдвиговой устойчивости на разных структурно-масштабных уровнях. В многоуровневой самоорганизации структурно-фазовых превращений принципиально важную роль играют все типы интерфейсов в нагруженном твердом теле, на которых возникает «шахматное» распределение растягивающих и сжимающих нормальных и касательных напряжений.
Дальнейшее развитие физической мезомеханики позволило адекватно описать возникновение «шахматного» распределения нормальных и касательных напряжений противоположных знаков на поверхностях и интерфейсах нагруженного твердого тела, распространение одиночных и двойных спиральных волн на поверхности, квазипериодическое растрескивание и отслаивание упрочняющих и термобарьерных покрытий и др.
Очень важным фундаментальным результатом физической мезомеханики явилась теоретически и экспериментально обоснованная концепция о том, что в основе возникновения нелинейных волн локализованного пластического течения лежит самоорганизация каналиро-ванных на мезомасштабном уровне потоков локальных структурных превращений, развивающихся на нано-масштабном уровне, и аккомодационных полос сброса на макромасштабном уровне [1, 2].
Эта концепция представляет методологическую основу для описания процессов самоорганизации локальных структурных превращений в деформируемом твердом теле на различных структурно-масштабных уровнях и построения многоуровневых моделей деформации и разрушения твердых тел.
Механика сплошных сред, хорошо зарекомендовавшая себя в решении прикладных задач теории упругости на макромасштабном уровне, изначально не ставила своей целью описание сред со структурой. Сложность применения методов механики сплошных сред к компьютерному моделированию динамики изменения внутренней структуры материала связана с рядом факторов. Одним из этих факторов является то, что в данных задачах должно рассматриваться не только наличие внутренней структуры, влияющей на перенос энергии, но и динамическая перестройка как самой структуры, так и конфигурации внутренних границ раздела.
Как правило, существующие на данный момент методы механики сплошных сред учитывают наличие
внутренней структуры путем задания так называемой «нелинейности свойств», т.е. пространственного разброса параметров моделируемой среды [11].
Однако при решении задач физической мезомехани-ки требуется не только явно учитывать внутреннюю структуру материала путем задания различных свойств отдельных элементов (как это традиционно делается в численных методах механики сплошных сред), но и учитывать иерархический характер взаимодействия между различными структурными элементами, а также трансформацию исходной структуры в процессе пластического формоизменения твердого тела.
Поскольку классическая механика сплошных сред не рассматривает внутреннюю структуру среды по определению, для имитации наличия структуры используются порой искусственные и зачастую крайне сложные математические приемы.
В силу невозможности явного задания специфики взаимодействия между отдельными элементами среды, в механике сплошных сред задаются закономерности изменения свойств элемента, как правило, полученные из экспериментальных данных о поведении материала на макроуровне. Так, например, закон изменения напряжения от деформации каждого локального элемента среды задается путем введения функции отклика, являющейся аналогом ст-е-диаграммы на макроуровне, а такие параметры, как предел текучести для модельного учета деформационного упрочнения, берутся из статистических закономерностей макроуровня (например, закон Холла-Петча).
Таким образом, для моделирования деформации материалов со сложной структурой, имеющей тенденцию к многоуровневой структурной трансформации, наряду с классическими методами сплошных сред необходимо привлекать класс методов, позволяющих в явном виде задавать как структуру моделируемой среды, так и законы взаимодействия между структурными элементами. Это так называемые дискретные подходы. Наиболее известными дискретными методами являются методы молекулярной динамики и клеточных автоматов.
Методы молекулярной динамики весьма эффективны при моделировании поведения систем на атомном уровне и, следовательно, жестко привязаны к потенциалу межатомного взаимодействия, нахождение вида которого является самостоятельной и трудоемкой задачей. Существенную сложность представляет моделирование реальных систем ввиду огромного числа атомов, поведение каждого из которых в рамках методов молекулярной динамики описывается по раздельности. Несмотря на множество интересных и значимых результатов, полученных с помощью методов молекулярной динамики [12, 13] и показавших хорошие перспективы применения этих методов, данные подходы нельзя считать достаточными для теоретического описания поведения материала как иерархически организованной среды.
Метод клеточных автоматов, изначально разработанный для моделирования передачи импульсов головного мозга, за последние полвека существенно расширил область своих применений. Он успешно используется для решения широкого спектра задач — от газодинамики до проблем роста костной ткани, взаимодействия клеток организма, кинетики распространения волн горения и образования дендритных структур при кристаллизации, перестройки поликристаллических структур в процессе холодной прокатки сталей и др. [14-19].
Клеточными автоматами принято называть сети из активных элементов, меняющих свое состояние в дискретные моменты времени по определенному набору правил в зависимости от того, какими были состояния самого элемента и его ближайших соседей по сети в предыдущий дискретный момент времени.
Благодаря «абстрактности» клеточного автомата, поведение которого регулируется исключительно логикой, появляется широкий диапазон возможностей применения этого подхода. Это связано с тем, что логические условия перехода элемента автомата из одного состояния в другое можно представить в виде произвольного отображения Rn ^ N, т.е. существует возможность конструировать правила переключения исходя из любого набора физических, химических или иных параметров, значения которых влияют на эволюцию системы.
Тем не менее, все вышеперечисленные подходы, так или иначе, оставались «одноуровневыми». Каждый масштабный уровень, равно как и каждый метод моделирования на этом уровне, характеризуется своими законами и своим языком описания. Наибольшую сложность представляет самосогласование этих методов, или возможность передачи информации с одного масштабного уровня на другой. Наиболее перспективным на данный момент метаязыком описания многоуровневых систем и процессов является физическая мезомеханика.
В связи с этим создание междисциплинарных многоуровневых подходов к описанию процессов распространения ротационных волн пластической деформации, связанных со структурными трансформациями на разных масштабных уровнях, является наиболее актуальной задачей в современной физике деформируемого твердого тела.
В настоящей работе проведено экспериментальное и теоретическое исследование явления фрагментации поверхностного слоя в условиях трибосопряжения. Предлагается гибридный метод моделирования поведения образца, подвергающегося внешнему механическому воздействию (растяжению, сжатию и др.). В рамках предлагаемого гибридного подхода образец представляется в виде активных элементов микромасштабного уровня. Этот подход позволяет моделировать трансформацию структуры на микроуровне: возникно-
вение очагов неупругой деформации, зародышей пор, трещин и т.д.
2. Фрагментация поверхностного слоя в трибо-контакте и формирование частиц износа
Поверхность и примыкающий к ней слой играют очень важную роль в деформации и разрушении твердых тел. Поэтому они выделяются в отдельную структурную подсистему, которой отводится особая роль в процессах диссипации подводимой энергии и релаксации напряжений. При трении и изнашивании роль поверхностных слоев значительно возрастает в связи с тем, что через них происходит передача механической энергии при контактном взаимодействии в объем материала и локализуются многие процессы трибологичес-кого характера. Геометрическая и временная дискретность фактической площади контакта определяет плотность и характер распределения источников темпера-турно-силовых возмущений в зонах фрикционного взаимодействия. В результате поверхностные слои испытывают различную степень и вид деформации в зависимости от уровня силовых и тепловых нагрузок, а также своего расположения относительно фрикционного контакта. Наибольшую деформацию испытывают поверхностные слои, непосредственно примыкающие к поверхности изнашивания. Площадь реального пятна контакта очень мала, что определяет высокие контактные напряжения в поверхностном слое. При относительном скольжении пар трения пятна перемещаются, существуя в течение небольшого промежутка времени. Температура в пятне контакта резко повышается, а в микрообъеме зоны пятна касания, находящемся под действием нормальных и сдвигающих напряжений, возникают циклические знакопеременные напряжения, сжимающие перед пятном и растягивающие за ним. Эти напряжения приводят к механическому и термическому усталостному разрушению поверхностного слоя.
Высокая степень пластической деформации в поверхностных слоях приводит к развитию полей напряжений, деформаций и скоростей смещений с очень высокими градиентными характеристиками. В поверхностном слое развивается вихревое пластическое течение, которое создает в материале высокие моментные напряжения. Они вызывают фрагментацию исходной структуры поверхностного слоя на мезомасштабном уровне и формирование границ раздела, параллельных плоскости трения. По мере развития процесса изнашивания происходят уменьшение размера фрагментов и формирование фрагментированной структуры. В верхних слоях, примыкающих к поверхности трения и деформирующихся в условиях циклического нагружения «растяжение - сжатие», процесс фрагментации протекает наиболее активно, а формирующаяся при этом фрагментированная структура принимает относительно
равноосную форму. Процесс фрагментации завершается формированием интерфейса между тонким сильно дефектным поверхностным слоем и основным материалом.
Исследование структуры поверхностного слоя стали 20Х18Н10Т после нагружения в парах трения методом электронной микроскопии высокого разрешения показало наличие муарового узора, возникающего из-за наложения изображений разориентированных в плоскости фольги микрообластей (муар вращения) (рис. 1, а, б). Анализ картин муара показывает, что разориентация конгломерата зерен А относительно смежного с ним материала может достигать ~16°-17°. Эти результаты экспериментально подтверждают формирование больше-угловых границ конгломератов зерен в плоскости фольги под действием локальных моментных напряжений. В условиях несовместности пластической деформации поверхностного слоя и основного материала на их интерфейсе возникает квазипериодическое «шахматное» распределение сжимающих и растягивающих нормальных напряжений. В условиях вихревого характера деформации при достижении критической величины растягивающих напряжений по границам фрагментов зарождаются микротрещины. Их рост и слияние приводят к отслаиванию фрагментов в виде частиц износа различного масштаба.
В физической мезомеханике распространение трещины рассматривается как поворотная мода деформации на макромасштабном уровне. В соответствии с законом сохранения момента количества движения, в материале, окружающем трещину, должны развиваться аккомодационные поворотные моды на более низких структурно-масштабных уровнях. В соответствии с этим появлению трещины в зоне ее зарождения должны предшествовать локализованная пластическая деформация, сопровождаемая материальными поворотами в поверхностном слое и вихревой аккомодационной упругой деформацией в подложке. Последняя обусловливает
фрагментацию материала как кристаллографические повороты на мезомасштабном уровне и формирование трансляционно-ротационных мезообъемов. Такие мезо-объемы должны локализованно формироваться впереди трещины и по ее берегам на протяжении всего процесса ее распространения. Естественно, скорость распространения трещины будет определяться скоростью формирования трансляционно-ротационных мезообъемов в объеме деформируемого материала. На скорость их формирования оказывают влияние напряженно-деформированное структурно-фазовое состояние материала. Изменяя эти факторы, можно управлять скоростью формирования трансляционно-ротационных мезообъемов и скоростью продвижения трещины. Таким образом, процесс изнашивания является вихревым и многомасштабным. Естественно, скорость процесса изнашивания должна зависеть от исходного состояния поверхностного слоя.
3. Описание метода возбудимых клеточных автоматов в рамках гибридного алгоритма
Численное моделирование на базе методов клеточных автоматов осуществляется согласно описанному ниже алгоритму:
1) формирование геометрии внутренней структуры образца;
2) задание распределения локальных свойств материала для каждого элемента автомата по всему объему (например, разброс значений прочности, упругих модулей, скорости деградации материала, плотности дислокаций и т.п.);
3) задание связей между элементами;
4) задание множества возможных состояний элементов и формулировка правил их переключения;
5) численная реализация перераспределения и динамики распространения возмущений между элементами в результате внешних воздействий.
Рис. 1. Формирование дефектной структуры (а) и муаровый узор (б) в тонком поверхностном слое стали 20Х18Н10Т (5 мкм от поверхности) в процессе изнашивания в паре трения (указаны значения углов разориентации в соответствующих областях), в — электронограмма к (а)
В рамках метода возбудимых клеточных автоматов (Stochastic Excitable Cellular Automata (SECA) method) моделируемый образец разбивается на элементарные объемы, каждый из которых моделируется с помощью клеточного автомата. Помимо состояния, элемент характеризуется параметрами, соответствующими моделируемому объему среды. Клеточные автоматы подразделяются на три основных типа: бистабильные, возбудимые, автоколебательные. Бистабильный автомат может находиться в одном из двух возможных состояний. Возбудимый автомат способен совершать последовательную цепочку переключений состояний под влиянием внешнего воздействия. Автоколебательные автоматы также обладают набором возможных состояний и пробегают их в отсутствие влияния извне. Исходя из особенностей моделируемых процессов, предполагающих распределение тепловой и механической энергии в материале, в качестве инструмента моделирования был выбран возбудимый клеточный автомат. Каждый такой автомат характеризуется определенным набором соседей на первой координационной сфере, а также числовыми параметрами, соответствующими материалу, содержащемуся в моделируемом объеме пространства, такими как модуль упругости, плотность, модуль сдвига, плотность дислокаций, удельная теплопроводность, удельная теплоемкость, коэффициент температурного расширения и др. При взаимодействии с соседними автоматами могут меняться тепловая и механическая составляющие энергии, а значит, и связанные с ними физические параметры (температура, энтропия, напряжение, деформация, плотность и т.д.).
В зависимости от специфики моделируемого процесса выбираются параметры автомата, изменение которых может привести к переключению состояния данного автомата. Например, при исследовании тепловых процессов такие состояния могут отвечать за агрегатные состояния вещества, характеризуемые определенными интервалами температуры и значениями теплоемкости, теплопроводности, плотности и т.д. При моделировании процессов, связанных с притоком и перераспределением механической энергии, целесообразно связать состояния автомата со стадиями нагружения твердого тела: упругая деформация, неупругая деформация, пред-разрушение, разрушение. Каждая из этих стадий может быть охарактеризована определенными значениями таких параметров автомата, как плотность, модуль упругости, модуль сдвига и др. Переключение между состояниями может осуществляться как детерминированно, так и стохастически. В зависимости от такого типа переключения автоматы делятся на детерминированные и вероятностные. В рамках предлагаемого SECA-метода клеточный автомат является вероятностным, т.е. каждому значению исследуемого параметра соответствует определенный набор значений вероятности переключения в соответствующее состояние.
Важным фактором при моделировании является выбор типа упаковки клеточных автоматов. Каждый из этих типов характеризуется собственным набором направлений распространения энергии и количеством соседей на каждой координационной сфере. Вообще говоря, тип упаковки автоматов может не совпадать с типом упаковки атомов в кристаллической решетке моделируемого материала. В рассматриваемой модели можно выбрать либо простую кубическую, либо гексагональную плотную упаковку автоматов. Данный выбор зависит от специфики моделируемого процесса. Так, например, при решении задачи о вдавливании инден-тора целесообразно выбрать гексагональную плотную упаковку, т.к. она позволяет более точно моделировать распространение полусферического фронта возмущений. В случае моделирования потоков тепловой и механической энергии, направленных от плоских равномерно нагружаемых границ вглубь образца, более удачным представляется выбор простой кубической упаковки, структура которой подходит для моделирования плоских деформационных и температурных фронтов. Варьирование этих двух типов упаковки автоматов полезно также при исследовании корректности модели того или иного физического процесса: достаточная близость результатов, полученных в разных упаковках, говорит о правильном построении модели.
В рамках предлагаемого подхода клеточный автомат обладает следующими параметрами, относящимися к материалу, содержащемуся в моделируемом объеме: модуль упругости, модуль сдвига, плотность, коэффициент термического расширения, удельная теплоемкость, удельная теплопроводность, совокупная энергия дефектов материала, максимальная скорость отклика среды на внешнее механическое воздействие, полная энергия, механическая энергия, тепловая энергия, температура.
Помимо этих параметров, как уже было сказано выше, важной характеристикой клеточного автомата является тип упаковки автоматов. В рамках предлагаемого метода реализованы возможности построения сети в простой кубической и гексагональной плотной упаковках. При кубической упаковке на первой координационной сфере автомата располагаются 6 соседей (элементарная ячейка представляет собой куб), при гексагональной — 12 (элементарная ячейка имеет форму додекаэдра). Так как при моделировании механического нагру-жения учитываются потоки энергии только от автоматов на первой сфере, подобное различие приводит к тому, что при гексагональной упаковке моделируемый материал обладает большим количеством направлений распространения энергии из локального мезообъема, чем при кубической упаковке автоматов. Однако при моделировании плоских фронтов механической энергии кубическая упаковка имеет преимущество перед гексагональной, т.к. позволяет повысить точность и быст-
родействие алгоритма вычисления потоков энергии в направлениях, параллельных и перпендикулярных фронту.
4. Моделирование переноса фронта неупругой деформации
Предлагается математическая модель деформации нагруженного твердого тела, в основе которой лежит положение о том, что процесс деформации, по сути, есть следствие перераспределения энергии между различными структурными элементами твердого тела и трансформации различных ее частей друг в друга [20, 21].
В рамках SECA-метода моделируется распределение упругой энергии в образце, подвергающемся внешнему механическому воздействию (растяжению, сжатию и др.). Моделируемый образец представляется в виде клеточного автомата — сети элементов. Данная сеть разделена на кластеры, каждый из которых отвечает за отдельное зерно, обладающее собственной ориентацией кристаллической решетки, характеризуемой углами Эйлера ф, п.
Входными параметрами модели являются начальные значения гидростатического давления, плотности и температуры каждого элемента {р0, р0, Т°, 0 <1<1- 1}, I — общее число элементов. Начальная упругая энергия г-го элемента вычисляется следующим образом:
,0 ^Ка
2У0
(1)
Энергия г-го элемента на каждом п-м временном шаге зависит от его энергии и энергии каждого к-го соседнего элемента на 1-й координационной сфере г-го автомата (0 < к<К - 1, К — число элементов на 1-й сфере) на (п - 1)-м временном шаге и вычисляется следующим образом:
А = д"—1 + ? ДА", (2)
к=0
где ДА"- — изменение упругой энергии г-го элемента в результате взаимодействия с к-м соседом на п-м шаге
по времени, которое рассчитывается следующим образом.
1. Вычисляется напряжение ст"- на границе г-го элемента и его к-го соседа на (п - 1)-м шаге по времени как разница давлений, действующих со стороны каждого из рассматриваемой пары элементов:
п—1
-Р" •
_;п—1 п— 1
СТк = Рк -
(3)
2. По формуле Торнбулла рассчитывается скорость границы под действием напряжения ст"—1:
=—т"кСТк •
(4)
Здесь тл — подвижность границы между г-м элементом и его к-м соседом:
т.к = (т0).к ехР
0
кВТ.к
1
(5)
В данном случае (т0) л — максимальное значение подвижности, которое зависит от вида материала, содержащегося в каждом из этих элементов:
Р—(У —■гк )2/т)
(т0) к =
Щ
(6)
где У., Ук — значения модуля упругости г-го и к-го элементов; с — эффективная скорость отклика среды на внешнее механическое воздействие; кв — постоянная Больцмана; Qik — энергия границы между г-м элементом и его к-м соседом на (п - 1)-м шаге по времени:
Qik = УHAGB'
9 — 9.-12! — 1п ^0
HAGB
11
HAGB
(7)
где Уддев — максимальная энергия границы, соответствующая максимальному углу разориентации кристаллической решетки; 9НА^ — максимальный угол разориентации; 9., 9к — эффективные углы ориентации решетки г-го и к-го элементов; Тл — значение температуры на рассматриваемой границе, зависящее от температур данных элементов Т и Тк следующим образом:
Т.к ="
Т + Т
(8)
Ориентация
кристаллической
решетки
Элементы клеточного автомата
Атомное ядро
Свободные электроны
Граница зерна
Граница зерен
Свободные электроны
Область кривизны
Поток механической энергии
Рис. 2. Схематические изображения зерна с кристаллической решеткой (а) и соседних зерен с разными ориентациями решетки (б)
Поликристаллический образец
Ориентация кристаллической решетки зерна
Рис. 3. Схема тройного стыка зерен в поликристаллической структуре
Иллюстрация разориентации кристаллических решеток между зернами приведена на рис. 2. На рис. 3 приведена схема тройного стыка зерен в моделируемой поликристаллической структуре.
В данном случае значение подвижности является величиной, обратной удельному (по объему) импульсу силы реакции материала, содержащегося в соседнем элементе.
Далее вычисляются объемная доля материала, переместившегося в соседний элемент Д^, и изменение механической энергии г'-го элемента в результате взаимодействия с к-м соседом ДАП:
ду;
vn-1дt
гк _ \ГП _
~уГ ~Кк ~~
ДАПк =дСПк ^У-
(9) (10)
Здесь Д — величина временного шага; 1С — размер элемента; УС — объем элемента.
Затем вычисляются полное относительное изменение объема г'-го элемента П и приток механической энергии ДА; г'-го элемента в результате взаимодействия со всеми К соседними элементами:
к
(11)
ДУп к дуп
ДСП =ДУ- =1ДУк =ЕДС;к,
ус к=1 УС к=1
ДАП = Е МП -
к=1
(12)
Здесь ДуП > 0, если объем г'-го элемента увеличился, в противном случае ДуП < 0. Значения относительного объема материала СП и полной механической энергии
А. П г'-го элемента на п-м временном шаге вычисляются следующим образом:
С П = С П-1+ДС П, (13)
АП = АП-1 +ДАП. (14)
Зная значение плотности рП-1 и изменение объема ДуП, можно вычислить изменение массы ДцП и ее новое значение цгП для г'-го элемента:
дцП =-рП-1 уП =-РП-1 УсДС П, (15)
цП = цП-1+ДцП - (16)
Новые значения плотности материала, содержащегося в г'-м элементе, рП и его гидростатического давления рП определяются выражениями
= цП =ь"-1+ДцП =рП-1+ДцП р; у у р; + у
р = — -
' у
(17)
(18)
Значения {ргП, 0 < г < I - 1} и (ДП, 0 < г < I - 1} становятся входными параметрами для (п + 1)-го временного шага. Выходными параметрами моделируемого образца являются конечные значения гидростатического давления элементов {р?, 0 < г < I - 1}, где N — число временных шагов.
Необходимо отметить, что в рамках излагаемого подхода граничные элементы делятся на три типа: «жесткие», «мягкие» и «промежуточные». Граничные элементы «жесткого» типа при взаимодействии с внутренними элементами образца не изменяют собственную энергию. Влияние таких элементов друг на друга не моделируется. «Мягкие» граничные элементы взаимодействуют со всеми соседями на первой координационной сфере, как внутренними, так и граничными, при этом их энергия меняется в соответствии с соотношениями для потока энергии. Наконец, элементы промежуточного типа не взаимодействуют с соседними граничными элементами, изменяя энергию в результате влияния внутренних соседних элементов на первой сфере.
5. Алгоритм вычисления локальных моментов сил и величины угла кручения при формировании кривизны в материале
Алгоритм расчета локальных моментов сил, изложенный в [22], был развит для случая наличия потоков механической энергии с учетом разворота и кручения материала в активном элементе клеточного автомата. Для случая гексагональной плотной упаковки элементов клеточного автомата реализован алгоритм распределения механической энергии по элементам, с помощью которого на каждом временном шаге для каждого г'-го элемента можно вычислить значение главного напряже-
ст,,, =
Рис. 4. Схема алгоритма вычисления главного напряжения и момента силы
ния ст. исходя из значений гидростатического давления данного элемента и его соседей по 1-й координационной сфере. Для этого вычисляются скалярные значения напряжений, действующих на г-й элемент со стороны каждого из его соседей (рис. 4):
ст = |Рк — р, 8(",к) =1 (19)
1к <0, 8(I, к) = 0. Здесь стк — гидростатическое давление к-го элемента на 1-й координационной сфере г-го элемента (для используемой гексагональной плотной упаковки число таких элементов равно 12); 8(г, к) — функция принадлежности к-го элемента моделируемой сети элементов. Значение этой функции равно 1, если элемент принадлежит рассматриваемому множеству, и равно 0, если не принадлежит. В нашем случае, если для всех к от 1 до 12 8(г, к)= 1, то г-й элемент является внутренним, в противном случае — граничным.
Затем вычисляется скалярное значение главного напряжения:
1 12
ст 1 = —Ест.к •
12 к=1
(20)
Значение вектора момента силы г-го элемента М. находится следующим образом.
1) Вычисляются координаты векторов, направленных от центра г-го элемента к центрам каждого из его соседей на 1-й координационной сфере (рис. 4):
Гк = гк - Г1, (21)
где г. — радиус-вектор центра г-го элемента; Гк — радиус-вектор центра к-го элемента на 1-й координационной сфере г-го элемента.
2) Для каждого к-го элемента, являющегося внутренним, определяются индексы I элементов, лежащих на 1-й координационной сфере г-го и к-го элементов.
3) Вычисляются скалярные значения напряжений, действующих со стороны каждого к-го элемента на каждый его 1-й элемент, лежащий внутри образца:
Рк — Р1 > 8(к>1) = 1
=0, 8(к, I) = 0.
4) Определяется вектор, направленный от центра каждого к-го элемента к центру каждого его 1-го внутреннего элемента:
Гк1 = Г1 - гк • (23)
5) Вычисляется вектор силы , действующей со стороны каждого к-го элемента на каждый его 1-й элемент, лежащий внутри образца:
• 5 • ст
1к,
(24)
где 51 — площадь поверхности контакта соседних элементов.
6) Определяется искомое значение вектора момента силы г-го элемента:
М. =ЕЕ[Гк, % ], (25)
к I
где [г.к, % ] — векторное произведение векторов, полученных по формулам (23) и (24).
Приращение локального момента силы влечет за собой формирование кривизны в материале. Кривизна, в
0
Рис. 5. Расчетная схема стержня круглого поперечного сечения при кручении (а) и вычисление изменения плотности дефектов (б)
свою очередь, инициирует потоки дефектов вдоль границ раздела на разных масштабных уровнях. Под действием вектора момента силы происходит закручивание материала в активном элементе на определенный угол. Этот материальный поворот приводит к искажению кристаллической решетки путем формирования множественных линий скольжения. Более того, определение угла кручения материала дает возможность количественно оценить изменение объемной доли дефектов. Для этого можно воспользоваться приближением модели кручения стержня (рис. 5, а).
Искомый угол у вычисляется следующим образом: ССХ - ББ1 г(ф + dф) - гф г^ф
У = -
БС
dx
dx ,
(26)
0 = ,
dx
где 0 — относительный угол закручивания.
Воспользуемся приведенной выше моделью. Угол поворота элемента у вычисляется следующим образом:
у = г0, (27)
где г — радиус элемента. Относительный угол закручивания 0 определяется по формуле 2 М
0
Gп г 4
(28)
где М = |М|; G — модуль сдвига. Таким образом, 2 М
Т=БП7' <29)
Изменение числа дефектов в элементе при его повороте на угол у определяется исходя из того, какое количество элементарных ячеек кристаллической решетки попало в сектор, соответствующий этому углу (рис. 5, б):
ДМ d = У),
^ се11
(30)
где (у) — площадь сектора, соответствующего углу у; £се11 — площадь сечения элементарной ячейки решетки кристаллической плоскостью, проходящей через центр ячейки; к — коэффициент пропорциональности, зависящий от вида упаковки атомов и материала. Из «первых принципов» известно, что каждый материал характеризуется своим значением энергии межатомного взаимодействия, а значит, коэффициент к для различных материалов оказывается различным. Но он зависит не только от значения данной энергии, но и от многих других факторов.
Величина площади сектора £.есДу) определяется выражением
$»<*( У)
у г
(31)
Рис. 7. Компоненты локальных моментов сил (а, б) и распределение деформаций (в) для образца с размером зерна 10 мкм
Рис. 8. Компоненты локальных моментов сил (а, б) и распределение деформаций (в) для образца с размером зерна 30 мкм
Подставляя выражения (29) и (31) в (30), получаем: кМ (32)
ANd =
GnrSc
cell
Изменение плотности дефектов в элементе Дра вычисляем как отношение к площади сечения элемента клеточного автомата
Д ДЖ- кМ
ДРа =-- =-. (33)
5са Gпr5ce115ca
Полученное соотношение позволяет вычислять изменение плотности дефектов в материале при формировании кривизны под действием локального момента силы.
6. Моделирование фрагментации поликристаллической структуры в поле локальных моментов сил
SECA-методом проведены численные эксперименты по механическому нагружению поликристаллических образцов с различными размерами зерна. Каждый образец моделировался посредством клеточного автомата с ГЦК-упаковкой активных элементов размером 1 мкм. Линейные размеры образцов составляли
100х 10х 100мкм. Механические и физические параметры моделируемого материала, содержащегося в активных элементах, соответствовали алюминию. В первом образце средний размер зерна составлял 10 мкм, во втором — 30 мкм. Общее время нагружения — 100 мкс, величина шага по времени — 1 нс. Начальные значения температуры всех элементов для каждого образца были равны 300 К, начальные значения деформаций и напряжений равнялись нулю. Структуры образцов и схема нагружения представлены на рис. 6.
Результаты численных экспериментов показали, что в процессе нагружения поликристаллических образцов вдоль межзеренных границ формируются области экстремальных значений локального момента силы, которые приводят к кристаллографическому повороту материала в теле конгломерата зерен. Вид таких областей приведен на рис. 7 и8.
Для детального исследования процесса кристаллографического поворота конгломерата зерен было построено распределение десятичных логарифмов моментов сил на поверхности нагруженного крупнокристаллического образца (рис. 9). Видно, что отдельные зерна в конгломерате фрагментируются областями со знаком
100 x, мкм
Рис. 9. Зеренная структура исходного образца (а) и фрагментация зерен в конгломерате (1, 2, 3) в распределении десятичных логарифмов компонент Xуглов разворота sign(Y;tух| + 6) на лицевой грани образца в момент времени 100 мкс (б)
поворота, противоположным знаку разворота конгломерата зерен.
7. Выводы
Развит метод возбудимых клеточных автоматов для моделирования формирования пространственной кривизны на межзеренных границах в поликристалле при возникновении неоднородного векторного поля локальных моментов сил. Метод существенно расширен посредством нового алгоритма расчета величин кривизны и кручения материала в поле локальных моментов сил. На основе данного метода проведено исследование процесса формирования градиентного поля моментов сил на границах зерен в поликристалле. Показано, что, наряду с модуляцией значений момента силы и формированием кривизны разного знака вдоль межзеренной границы, под действием такого градиентного механического поля отдельные конгломераты самосогласованно деформирующихся зерен испытывают кристаллографические повороты. Это сопровождается фрагментацией отдельных зерен конгломерата с поворотом обратного знака. Результаты численного эксперимента хорошо согласуются с экспериментально наблюдаемым вихревым движением конгломератов зерен в поверхностных слоях трибоконтакта. Предлагается механизм формирования частиц износа на поверхностях трения.
Работа выполнена при финансовой поддержке проекта 111.20.1.1 СО РАН и грантов Президента РФ для поддержки ведущих научных школ (НШ-6116.2012.1) и РФФИ (проект № 11-01-00646).
Литература
1. Панин В.Е., Егорушкин В.Е., Панин A.B. Эффект каналирования пластических сдвигов и нелинейные волны локализованной пластической деформации и разрушения // Физ. мезомех. - 2010. -Т. 13. - № 5. - С. 7-26.
2. Панин В.Е., Егорушкин В.Е., Панин A.B. Физическая мезомеханика
деформируемого твердого тела как многоуровневой системы. I. Физические основы многоуровневого подхода // Физ. мезомех. -2006. - Т. 9. - № 3. - C. 9-22.
3. Егорушкин В.Е. Динамика пластической деформации. Волны локализованной пластической деформации в твердых телах // Физическая мезомеханика и компьютерное конструирование материалов / Под ред. В.Е. Панина. - Новосибирск: Наука, 1995. -T.1.- C. 50-77.
4. Егорушкин В.Е. Динамика пластической деформации. Волны лока-
лизованной пластической деформации в твердых телах // Изв. вузов. Физика. - 1992. - № 4. - С. 19-41.
5. Zuev L.B., Barannikova S.A. Evidence for the existence of localized plastic flow auto-waves generated in deforming metals // Natural Science. - 2010. - V. 2.- No. 5. - P. 476-483.
6. Зуев Л.Б., Данилов В.И., Горбатенко B.B. Автоволны локализован-
ной пластической деформации // ЖТФ. - 1995. - Т. 65. - № 5.-С. 91-103.
7. Zuev L.B. Wave phenomena in low-rate plastic flow in solids // Ann. Phys. - 2001. - V. 10. - No. 11-12. - P. 965-984.
8. Панин B.E., Егорушкин B.E., Панин A.B. Нелинейные волновые процессы в деформируемом твердом теле как многоуровневой иерархически организованной системе // Успехи физических наук. - 2012. - Т. 182. - № 12. - С. 1351-1357.
9. Панин A.B. Нелинейные волны локализованного пластического те-
чения в наноструктурных поверхностных слоях твердых тел и тонких пленках // Физ. мезомех. - 2005. - Т.8.- № 3. - С. 5-17.
10. Физическая мезомеханика и компьютерное конструирование материалов / Под ред. В.Е. Панина. - Новосибирск: Наука, 1995. -Т. 1. - 298 с., Т. 2. - 320 с.
11. Романова B.A., Балохонов P.P. О роли внутренних границ раздела в процессах формирования мезоскопического деформационного рельефа на свободной поверхности нагруженных металлов // Физ. мезомех. - 2010. - Т. 13. - № 4. - С. 35-44.
12. Abraham F.F. Lecture Notes in Physics. - 2001. - V. 567. - 436 p.
13. Krivtsov A.M. Molecular dynamics simulation of impact fracture in polycrystalline materials // Meccanica. - 2003. - V. 38. - No. 1. -P. 61-70.
14. Лоскутов A.Ю., Михайлов A.C. Введение в синергетику. - M.: Наука, 1990. - 272 с.
15. Wolfram S. Cellular automata as models of complexity // Nature. -1984. - V. 311. - P. 419-424.
16. Полак Л.С., Михайлов A.C. Самоорганизация в неравновесных физико-химических системах. - M.: Наука, 1983. - 286 с.
17. Bmep П., Розенблют A. Проведение импульсов в сердечной мышце // Кибернетический сборник. - Вып. 3. - M.: ИЛ, 1961. - С. 356.
18. Kroc J. Application of сellular automata simulations to modelling of dynamic recrystallization // Lecture Notes in Computer Science. -2002. - V. 2329. - P. 773-782.
19. Godara A., Raabe D. Mesoscale simulation of the kinetics and topology of spherulite growth during crystallization of isotactic polypropylen (iPP) by using a cellular automaton // Modelling Simulation Materials Science Engineering. - 2005. - V. 13. - P. 733-751.
20. Моисеенко Д.Д., Максимов n.B. Распределение напряжений и деформаций на интерфейсе «поверхностный слой - подложка»: моделирование на основе стохастического подхода // Физ. мезомех. -2005. - Т. 8.- № 6. - С. 89-96.
21. ПанинB.E., Моисеенко Д.Д., Максимов n.B., Панин A.B. Физическая мезомеханика деформируемого твердого тела как многоуровневой системы. III. Неупругий предвестник зарождения пластического сдвига // Физ. мезомех. - 2006. - Т. 9. - № 5. - С. 5-15.
22. Панин B.E., Сергеев B.П., Моисеенко Д.Д., Почивалов Ю.П. Научные основы формирования теплозащитных и износостойких многослойных покрытий системы Si-Al-N/Zr-Y-O // Физ. мезомех. -2011. - Т. 14. - № 6. - С. 5-14.
Поступила в редакцию 24.09.2012 г.
Сведения об авторах
Моисеенко Дмитрий Давидович, к.ф.-м.н., снс ИФПМ СО РАН, mdd@ispms.tsc.ru Почивалов Юрий Иванович, к.ф-м.н., внс ИФПМ СО РАН, pochiv@ispms.tsc.ru Максимов Павел Васильевич, к.ф.-м.н., нс ИФПМ СО РАН, mpv@ispms.tsc.ru
Панин Виктор Евгеньевич, д.ф.-м.н., академик РАН, научн. рук. ИФПМ СО РАН, paninve@ispms.tsc.ru
