CC BY
37
УДК 539.3
Прямая физическая уируговязкоиластическая модель: приложение к исследованию деформирования монокристаллов
П.В. Трусов1, А.Ю. Янц1, Л.А. Теплякова2
1 Пермский национальный исследовательский политехнический университет, Пермь, 614900, Россия 2 Томский государственный архитектурно-строительный университет, Томск, 634003, Россия
Представлены постановка и результаты моделирования одноосной осадки монокристалла алюминия технической чистоты, приведено сравнение с экспериментальными данными. Для решения задачи использована прямая модель первого типа, основанная на методе конечных элементов. Поведение материала описывалось упруговязкопластической моделью физической теории пластичности, которая позволяет явно учитывать сдвиги по кристаллографическим плоскостям. Основной особенностью данной работы является физически обоснованное описание геометрической нелинейности, связанной с вращением кристаллической решетки. Результаты моделирования показали, что исходный однородный монокристалл претерпевает фрагментацию на объемы, отличающиеся друг от друга различной интенсивностью скорости пластических сдвигов и ориентацией кристаллической решетки, что находится в удовлетворительном соответствии с данными экспериментальных исследований.
Ключевые слова: монокристалл, геометрическая нелинейность, физические теории пластичности, краевая задача, прямая модель, фрагментация, кривизна решетки, алюминий
DOI 10.24411/1683-805X-2018-12004
Direct elastoviscoplastic model: an application to the study of single crystal deformation
P.V. Trusov1, A.Yu. Yanz1, and L.A. Teplyakova2
1 Perm National Research Polytechnic University, Perm, 614990, Russia 2 Tomsk State University of Architecture and Civil Engineering, Tomsk, 634003, Russia
The uniaxial compression of a single crystal of commercial purity aluminum has been modeled, and the modeling results have been compared with experimental data. The problem is solved using the first-type direct model based on the finite element method. The material behavior is described by an elastoviscoplastic model of plasticity theory which explicitly accounts for shearing along crystallographic planes. The main feature of this study is a physically justified description of the geometric nonlinearity associated with crystal lattice rotation. The modeling results show that the original homogeneous single crystal is divided into volumes with different plastic shear rate intensity and lattice orientation, which is in satisfactory agreement with experimental data.
Keywords: single crystal, geometric nonlinearity, plasticity theories, boundary value problem, direct model, fragmentation, lattice curvature, aluminum
1. Введение
Устойчивый интерес к экспериментальным исследованиям поведения монокристаллов металлов и сплавов при различных термомеханических воздействиях сохраняется в течение нескольких десятилетий [1-6]. Данное обстоятельство обусловлено тем, что монокристалл является объектом с относительно простой внутренней структурой, что позволяет детально анализировать основные (в первую очередь дислокационные) механизмы неупругого деформирования. Наряду с экспери-
© Трусов П.В., Янц А.Ю., Теплякова Л.А., 2018
ментальными, в последние 15-20 лет интенсивно развиваются теоретические исследования, основанные на применении многоуровневых физических моделей деформирования монокристаллов [7-17]. Указанные модели деформирования монокристаллов позволяют выявить процессы и механизмы, реализующие неупругое деформирование, которые недоступны при экспериментальном исследовании, особенно внутри объема исследуемых образцов. Объединение усилий физиков и механиков, специализирующихся на разработке
физически обоснованных теорий, способствует совершенствованию как моделей для описания поведения монокристаллов, так и методов экспериментального исследования, что в итоге позволит уточнять и дополнять модели неупругого деформирования монокристаллов и приведет к возможности более точного описания поведения поликристаллов. Следует также отметить, что в ряде отраслей промышленности (в первую очередь, авиадвигателестроении) в последние 1015 лет предпринимаются значительные усилия по созданию монокристаллических деталей и конструкций, в связи с чем развитие теоретических и экспериментальных методов анализа поведения монокристаллов при различных термомеханических воздействиях востребовано и для исследования процессов изготовления и эксплуатации реальных монокристаллических конструкций.
Известно, что поведение монокристаллов при деформировании, эволюция дефектной структуры существенным образом зависят от взаимной ориентации характерных осей нагружения и кристаллической решетки, что связано с реализацией движения краевых дислокаций по различным наборам систем скольжения при отличающихся ориентациях. При этом общее количество таких систем зависит от типа кристаллита (например, для ГЦК-кристаллов — 12 систем скольжения, для ОЦК — 48) и полностью определяется его кристаллографией [18, 19]. В свою очередь, реализация скольжения по различным наборам систем скольжения приводит к существенным отличиям в эволюции напряженного состояния и дефектной структуры монокристалла.
При экспериментальном исследовании процессов нагружения монокристаллов выявлен ряд общих закономерностей. В работах [1, 20-22] показана существенная зависимость характера развития сдвиговой деформации от ориентации кристаллита относительно оси сжатия. В работах [20, 23-25] указывается на возникновение локализации сдвига в монокристаллических образцах. В [24] представлены результаты исследования процессов фрагментации монокристаллических образцов с ГЦК-решеткой при сжатии, приводящих к выделению в объеме изначально однородного монокристалла специфических зон, в которых активны различные наборы систем скольжения. Расположение и конфигурация указанных зон существенно зависят от краевых условий, таких как расположение и ориентация (относительно характерных направлений нагружения) свободных и контактных поверхностей монокристалла, условий на поверхностях контакта кристалла и штампа, взаимной ориентации граней кристаллического образца и кристаллографических систем скольжения.
В связи со значительным влиянием на неоднородность деформирования граничных (и особенно контактных) условий формулировке последних уделена боль-
шая часть раздела 5. Контактные граничные условия формулируются как смешанные: кинематические для нормальной компоненты и силовые — для касательных составляющих. При этом предусмотрена возможность прилипания частиц контактной поверхности к плите пресса, возможность нарушения контакта и возможность скольжения (закон Зибеля-Амонтона-Ку-лона). В то же время в настоящей работе не ставилась цель специального исследования влияния контактных условий на поведение монокристаллического образца, эта задача составляет предмет отдельного рассмотрения.
Взаимные ориентации систем скольжения и свободных граней оказывают существенное влияние на величину и распределение пластических сдвигов в объеме кристалла. Наряду с этим, в экспериментах с увеличением степени деформации наблюдается искривление первоначально прямых следов полос выхода краевых дислокаций на поверхность, а также относительный поворот октаэдрических площадок {111}. Во всех указанных работах отмечается также, что при пластическом деформировании монокристаллов обнаруживается существенная неоднородность сдвиговых деформаций, а следовательно, и напряженно-деформированного состояния.
Известно, что пластическое деформирование кристаллов реализуется преимущественно путем движения дислокаций по кристаллографическим плотноупакован-ным плоскостям, при этом основной вклад в неупругое деформирование вносит движение краевых дислокаций [18, 19, 26-28]. В связи с этим процессы фрагментации объема монокристалла существенным образом зависят от кристаллогеометрии деформируемого образца и краевых условий. Под кристаллогеометрией монокристаллического образца будем понимать совокупность геометрических параметров, описывающих взаимные ориентации оси (или осей) нагружения, боковых граней и систем скольжения. Для различающихся взаимоориен-таций осей нагружения и кристаллической решетки реализуется скольжение по различным наборам систем.
Как было отмечено ранее, при деформировании монокристаллов имеет место неоднородность пластических сдвигов, дислокационных субструктур, напряжений и разориентаций одних объемов кристалла относительно других. В конечном итоге при деформациях порядка 3 % и выше в изначально однородных монокристаллических образцах возникает разделение на зоны, в которых эволюция внутренней структуры происходит по различным сценариям. В связи с этим ставится задача исследования процессов возникновения и эволюции локальных областей с различным состоянием внутренней структуры монокристаллического образца при неупругом деформировании для различных ориентаций осей нагружения, решетки и граней монокристалла.
Отметим, что в настоящей статье особое внимание уделено описанию процессов разориентации отдельных областей монокристалла при одноосной осадке, в связи с чем отдельно будет рассмотрена задача выделения квазитвердого движения из полного, которая имеет большое значение при рассмотрении геометрически нелинейных проблем [29, 30].
В настоящее время для описания поведения моно-и поликристаллических металлов и сплавов широко используется иерархический подход, основанный на рассмотрении материала как многоуровневой системы с присущими каждому уровню особенностями [31-33], при этом масштабы описания могут достигать микро-и наноразмеров. К моделям данного класса относятся также модели, основанные на физических теориях пластичности (упруговязкопластичности), которые позволяют явным образом включить в рассмотрение движение дислокаций по кристаллографическим системам скольжения, анизотропию упругих и неупругих свойств, которыми обладают монокристаллы. В настоящей работе для описания неупругого поведения отдельных кристаллитов используется упруговязкопластическая модель указанного класса [34], которая обладает рядом достоинств, например отсутствием неоднозначности выбора активных систем скольжения (которая присуща упругопластическим моделям) и относительной устойчивостью численного счета. Следует отметить, что при построении определяющих соотношений для описания процессов неупругого деформирования монокристаллов необходимо учитывать наличие геометрической нелинейности [33], включая учет поворотов кристаллической решетки при деформировании.
2. Модель для описания деформирования кристаллита
Остановимся вкратце на модели для описания деформирования кристаллита, основанной на физической теории упруговязкопластичности [34]. В качестве меры скорости полных деформаций используется градиент относительной скорости перемещения z (относительно подвижной системы координат, связанной с решеткой кристаллита), являющийся несимметричным и не зависящим от выбора системы отсчета (индифферентным) [35] тензором 2-го ранга [36]:
об
z = VvT =VvT-ю,
(1)
где V — оператор Гамильтона (набла-оператор), определенный в актуальной конфигурации; V — скорость перемещения относительно системы координат, вращающейся с угловой скоростью, определяемой тензором спина ю; V — скорость перемещения; ю — спин квазитвердого вращения жесткой системы координат, связанной с кристаллической решеткой материала, определенный ниже. Неголономная мера полных де-
формаций е определяется коротационным интегрированием введенной меры скорости деформации:
есг = е + е • ю - ю • е = Уут - ю, (2)
где (•)" — коротационная (индифферентная) производная. Более подробное рассмотрение данной меры, в том числе ее физический смысл, представлено в работе [37].
Для введенной меры скорости деформации принимается аддитивное разложение на упругую и неупругую составляющие:
г = ге + г1п. (3)
В данной модели полагается, что в начальный момент времени в объеме материала содержится достаточное количество краевых дислокаций, необходимое для начала процесса неупругого деформирования. Стоит отметить, что в модели не рассматриваются вопросы зарождения, движения и аннигиляции отдельных дислокаций; движение дислокаций описывается путем введения сдвигов по кристаллографическим системам скольжения [18, 19]. Каждая система скольжения характеризуется единичным вектором нормали плоскости скольжения п и направлением — нормированным вектором Бюргерса Ь, которые имеют фиксированные направления относительно кристаллографической системы координат на протяжении всего процесса. В каждый момент времени неупругая составляющая градиента относительной скорости перемещений определяется следующим образом:
= Е у k=i
(к )ь( к )n(k)
(4)
& (к)
где у4 ' — скорость сдвига по системе скольжения с номером Ь Скорость сдвига в системе скольжения определяется вязкопластическим соотношением [38]
(то) У/т _
Н(т(у) - т^)), у = 1, К, (5)
у (J) = У о
T(j)
где т( у), тСУ) — действующие и критические касательные напряжения в системе скольжения с номером j^; у 0 — скорость сдвига в системе скольжения при достижении действующим касательным напряжением значения критического напряжения; т — параметр скоростной чувствительности материала; К — общее количество систем скольжения в кристалле. Для нахождения критических касательных напряжений в системе скольжения с номером k может быть использовано эволюционное уравнение вида
J
& (к) _
= V f (к) = ^ f(j) j=1
(у'
(i)
у(i) ;а(Jp]), i, к = 1, N,
(6)
в котором каждое слагаемое ответственно за тот или иной механизм упрочнения. Стоит обратить внимание на то, что в настоящей работе упрочнение не рассматривается, т.е. для всех систем скольжения Tj,j) = const. Данное допущение принято с целью концентрации вни-
мания на анализе разориентации локальных объемов исходного однородного монокристалла.
В качестве определяющего соотношения принят анизотропный закон Гука в скоростной релаксационной форме:
осг = о + о • ю - ю • о = п-/т), (7)
где п — тензор четвертого ранга упругих свойств кристаллита.
3. Геометрическая нелинейность
В настоящее время одним из нерешенных вопросов нелинейной механики деформируемого твердого тела при рассмотрении процессов деформирования с большими градиентами перемещений является проблема выделения квазитвердого движения из полного, или разложение движения на деформационную и квазитвердую составляющие. Данный вопрос тесно связан с выбором коротационной производной, введение которой в определяющих соотношениях скоростного типа необходимо для удовлетворения принципа независимости от выбора системы отсчета. Следует отметить неоднозначность выбора коротационной производной, поскольку с формальной точки зрения существует множество мощности континуум способов разложения движения. В связи с этим главное отличие используемых разложений заключается в степени физической корректности. В частности, при построении моделей, основанных на физических теориях пластичности, используются два основных подхода к разложению движения. Первый основан на полярном разложении градиента места, второй — на разложении градиента скорости перемещений на симметричную и вихревую составляющие. Основным недостатком данных подходов [39] является отсутствие связи жесткой подвижной системы координат, отвечающей за квазитвердое движение, с выделенными на протяжении всего процесса кристаллографическими направлениями, характеризующими, в свою очередь, симмет-рийные свойства материала. В связи с этим указанные модели не имеют достаточно убедительного физического обоснования для их применения в описании вращения кристаллической решетки. В работах [30, 33, 39] показано также, что корректные определяющие соотношения для монокристаллов (т.е. анизотропных материалов) можно построить, если квазитвердое движение связано с решеткой, при этом компоненты тензора упругих характеристик в базисе, связанном с кристаллографической системой координат, остаются неизменными при любых наложенных жестких движениях. В данной работе используется подход [33, 39], основанный на «привязке» на протяжении всего процесса введенной жесткой подвижной системы координат, отвечающей за квазитвердое движение, к одному кристаллографическому направлению и содержащей его кристаллографической плоскости монокристалла.
Рассмотрим предлагаемый способ разложения движения более подробно. Введем лагранжеву кристаллографическую систему координат Оугу2у3 с базисом «вмороженную» в решетку. В отсчетной конфигурации примем эту систему декартовой ортогональной, совместив ее оси, например, для кубической неискаженной решетки с кристаллографическими направлениями [100], [010], [001]. При произвольном движении деформируемого кристалла этот триэдр будет претерпевать как повороты, так и искажения, т.е. его ортонор-мированность будет нарушаться. Заметим, что для кристаллических металлов и сплавов собственно искажения будут малыми.
Наряду с лагранжевой кристаллографической системой координат введем жесткую декартову ортогональную подвижную систему координат ОхХх3 с орто-нормированным базисом Ц., в отсчетной конфигурации совпадающую с кристаллографической системой координат и связанную с кристаллографической системой координат в течение всего процесса деформирования. Связь осуществим следующим образом: оси Оу1 и Ох1 будем считать совпадающими в каждый момент времени (вектор к1 направлен вдоль вектора q1); вектор к2 в каждый момент деформирования будем располагать в плоскости Оу 1у2 ортогонально вектору к1. Зная в каждый момент времени положение векторов к1, к2, легко определяется положение третьего базисного вектора подвижной системы координат: к3 = к1 х к2. Поскольку оси подвижной системы координ ат в данном случае связаны с кристаллитом, любое жесткое движение последнего будет воспроизводиться движением подвижной системы координат. Полагается, что пластические деформации не изменяют ориентации решетки, повороты и искажения последней определяются только упругой составляющей Ууе градиента скорости перемещений. Опуская выкладки, спин подвижной системы координат ю относительно условно неподвижной лабораторной системы координат будет выражаться соотношением
ю = кгк; = -(к2 VVе • М^к2-
- (к3 • Ууе • к! )Цк3 + (к2 -VVе • к! )к2^ -
- (к3 • Ууе •к2)к2к3 + (к3 • VVе • к1 )к3к1 +
+ (к3 •Ууе • к2)к3к2. (8)
Полученное выражение позволяет определить тензор спина квазитвердого движения кристаллита по известному значению упругой составляющей градиента скорости перемещений. При использовании данного подхода квазитвердое вращение кристаллической решетки будет определяться эволюцией пластических сдвигов в материале и градиентом полных скоростей перемещений.
4. Алгоритм вычисления кривизны решетки
Для определения разориентации объемов монокристалла в процессе деформирования необходимо вычисление значений кривизны и кручения кристаллической решетки. В каждой точке кристаллической решетки материала можно использовать кристаллографический лагранжев базис с координатами {у1}, введенный выше для определения спина решетки, векторы которого совпадают с выделенными осями в кристалле. При произвольном деформировании кристаллическая решетка испытывает вращение как жесткое целое и упругое искажение. Пластические деформации не приводят к искажению решетки. При этом упругие искажения решетки являются малыми по сравнению с неупругими деформациями, которым подвергается монокристалл. В связи с этим при определении кривизны будем пренебрегать незначительными упругими искажениями решетки, а принимать во внимание только вращение как жесткого целого. Для описания чисто вращательного движения используется модель ротации, описанная в предыдущем разделе, в которой вводится ортонормиро-ванная подвижная система координат с координатами {х'}, базисные векторы которой {к.} жестко связаны с одним кристаллографическими направлением и плоскостью. Ориентация подвижной системы координат относительно лабораторной системы координат с базисом {е.} и координатами {X'} определяется ортогональным тензором:
О = к 1е1,
0тп = ет ' кгег ' еп = ет 'кп = (9)
= кт 'кгег ' кп = ет ' кп-
В этом случае характеристики кривизны и кручения кристаллической решетки будут полностью определяться различием ориентаций подвижной системы координат в бесконечно близких точках кристаллической решетки. Следуя [40], тензор кривизны определяется соотношением
к = —2е :(От -(ОУ)),
(10)
где е — абсолютно кососимметричный тензор третьего ранга Леви-Чивиты. Выражение для компонент тензора кривизны кристаллической решетки в базисе лабораторной системы координат:
2 е: (от-(оу))=
к = —
-1 е 2
Ое jei
Отпетеп
дХ'
к ек
О,
дОп
дХ
к еУепек
= —1 е
е е е
'рдг р д г
о,
дОп
дХ
к еjenek
= —Iе о эа 2 дХ
к ерек,
к =— I е О ^Отп.
2 п т дхj ' Используя независимость векторов базиса лабораторной системы координат {ег} от координат {X1} и (9)2, преобразуем последнее соотношение:
к.. = —1 е- О дОтп=—1е. О д(ет •кп) =
у ~ та^та у ~ та^та
2 дХ 2 дXJ
= —1е О е дкп =—1е к -ее дк" (12)
~ та^та т у л та а т т ' (12)
2 дХ^ 2 дXJ
к =—1е к .^кп
у 2 тч 4 дхУ'
Тензор к в некоторой точке г кристаллической решетки описывает кривизну и кручение трех материальных решеточных линий Г, которые строятся как огибающие касательных к векторам Ц базиса подвижной системы координат в окрестности точки г. Диагональные компоненты к у тензора кривизны характеризуют кручение кривой 1', т.е. «градиент» (точнее, производную по направлению) угла поворота триэдра подвижной системы координат вокруг оси Ц вдоль этой же оси. Недиагональные компоненты к. характеризуют кривизну кривой 1' в плоскости, ортогональной вектору Ц. Иначе говоря, недиагональные компоненты к. (' ^ у) тензора к характеризуют изгибы кристаллографических осей решетки.
Для точного определения кривизны и кручения решетки необходимо непрерывное поле ориентаций подвижной системы координат, связанных с кристаллической решеткой. Однако при решении реальных задач известны ориентации подвижной системы координат в конечном множестве точек, в связи с чем вычисление кривизны возможно лишь приближенно. Основной сложностью при вычислении компонент тензора кривизны является вычисление компонент градиента орто-тонального тензора дОтп/ дХу.
В настоящей работе для решения краевой задачи используется метод конечных элементов. В связи с этим предлагается, используя некоторые функции формы, использовать конечно-элементную интерполяцию значений ортогонального тензора внутри каждого конечного элемента:
Оу (Г) = £ ф(р) (Г )О.) (г(р)), (13)
р
где ф( к) — функции формы конечного элемента; Оур) — значение компонент ортогонального тензора в р-м узле конечного элемента, при этом их определение является отдельной задачей, решение которой предложено ниже; г — радиус-вектор, лежащий в пределах конечного элемента; г(р) — радиус-вектор р-го узла конечного (11) элемента. Тогда значение компонент градиента Оу внутри конечного элемента вычисляется следующим образом:
dOm
dXj
(Г ):
Е ф(к) (Г)0% (Г(к)) =
= Е 0( к)
к
dXj к эф( к) (г) dXj
(к ).„(к К_
(14)
Таким образом, значение компонент тензора кривизны к внутри конечного элемента полностью определяется значением компонент ортогональных тензоров в узлах конечного элемента и функциями формы данного конечного элемента:
к j(Г)
к j(Г)
.-2
(Г)
(Е 0(к) дФ(к) (Г) ^ 0
dXJ
-1 €,.,
(
Е ф( p) МОЙ?
Е О
к
(к) эф(к) (г) ^ dXj
(15)
где использованы соотношения (11)-(14), г принадлежит конечному элементу.
Следует отметить, что при применении метода конечных элементов скорость вращения подвижной системы координат (определяемая моделью (8)), связанной с кристаллической решеткой, вычисляется в точках интегрирования, которые жестко связаны с материалом (т.е. имеют постоянные лагранжевы координаты) и находятся внутри конечного элемента. Отметим также, что при работе с ориентациями трехмерных объектов удобно использовать кватернионы, компоненты которых взаимно-однозначно связаны с компонентами ортогонального тензора [41]. Кватернион имеет четыре компоненты, которые выражаются через угол вращения ф и вектор и, направленный вдоль оси вращения:
I Ф . Ф • Ф • Ф q,- = < cos—, sin—u,, sin—u2, sin—u3 1 | 2 2 2 2 3
I Ф • Ф q = {cos—, sin—u
I 2 2
(16)
В процессе решения значения компонент кватерниона q вычисляются только в точках интегрирования, после чего (на каждом шаге решения) отдельно вычисляются значения компонент кватернионов в узлах конечного элемента сетки путем усреднения по смежным с узлом конечным элементам следующим способом. Для каждого узла N определяется массив точек интегрирования Iк тех конечных элементов, в которые входит рассматриваемый узел. Далее для этих точек интегрирования определяются весовые коэффициенты следующим образом:
s-1 ^ 1 Щ =~т' 5 = ¿7'
¡к i Ч
где ¡к — расстояние от рассматриваемого узла до &-й точки интегрирования (к е 1к). После этого вычисляется среднее значение ориентации кристаллической
решетки по весовым коэффициентам и значениям в смежных точках интегрирования:
q n = Е wq к , (17)
к
где q N, q к — кватернионы, определяющие ориентацию в узле с номером N и в смежных точках интегрирования 1к соответственно. Соотношение (17) является обобщением линейной интерполяции между двумя ориента-циями [41]. Для оценки точности данного способа получения узловых значений ориентации подвижной системы координат вычислялся угол несоответствия 8 между ориентацией qj в точке интегрирования с номером I, получаемой напрямую из модели материала, и ориентацией в точке интегрирования qj, найденной интерполяцией по узловым значениям:
¿1 j =Е ф^ )(ri) qN , (18)
N У '
где q N — кватернион, определяющий ориентацию в узле с номером N, входящем в конечный элемент с точкой интегрирования I, вычисленный по (17); r — радиус-вектор положения точки интегрирования с номером I. Известно, что последовательность двух поворотов, заданных кватернионами q,, q2, определяет поворот q = q2 ° q,, где ° — кватернионное умножение, результат которого также кватернион. В нашем случае можно вести кватернион q5, который совмещает ориентации qI и qI: qI = q5 ° qI. Тогда угол 8, на который задает вращение кватернион q5 (16) и есть искомое несоответствие. В численном эксперименте определялась величина
Д= max (5,-/ф-), i$t >0.001
где i — номер точки интегрирования; ф.—угол поворота подвижной системы координат кристаллической решетки в точке интегрирования i. При этом поиск максимума производился по всем точкам интегрирования, в которых поворот подвижной системы координат был не меньше 10-3 рад. В результате получено, что величина Д не превышает значение 15 %.
5. Прямая модель. Применение метода конечных элементов
Рассматривается задача одноосного сжатия монокристалла Al ориентировки D1 [24] (рис. 1), боковые грани совпадают с кристаллографическими плоскостями (001) и (110), а ось сжатия совпадает с кристаллографическим направлением (110). Монокристалл алюминия технической чистоты с размерами 6x3x3 мм зажат между двумя пуансонами P, и P2, между пунсонами и образцом имеет место трение, верхний пуансон движется поступательно со скоростью v = {0, 0, -V}, нижний пуансон неподвижен.
В момент времени t монокристаллическое тело занимает область О,' с границей S'. Необходимо опре-
x
Рис. 1. Схема расчетной области монокристалла ориентировки D1 (размеры указаны в мм) и координатные оси лабораторной системы координат, Р1, Р2 — верхний и нижний пуансоны, 51 = A1B1C1D1' S2 = ABCD, S3-6 — боковые поверхности
делить поля скоростей перемещений V и напряжений о в любой момент времени, удовлетворяющие уравнению равновесия (в скоростной форме) [29]:
V• о-V-(у-Уо) = 0, ге &, (19)
V — набла-оператор в текущей конфигурации, и определяющим соотношениям модели материала (1)-(8). Граница монокристалла в общем случае представи-ма в виде
= 5 &ее(') и 5соп*('), (20)
5 &ее(') — свободная граница; 5 соп*(') — граница контакта (смешанные граничные условия — силовые и кинематические). Поверхности ЛБСВ и Л1Б1С1В1 являются областями возможного контакта, на них задаются смешанные граничные условия. Следует отметить, что контактная граница 5соп1(') априори неизвестна. Для ее определения используется подход, подробно
описанный в [29], при этом сама граница контакта раз*
бивается на три участка: зона прилипания 5са, зона
* *
отлипания 5 ^ и зона скольжения 5 с8.
Ниже представлена полная совокупность граничных условий рассматриваемой краевой задачи деформирования монокристаллического образца, все компоненты определены в базисе неподвижной лабораторной системы координат:
1) верхний пуансон совершает поступательное прямолинейное движение с заданной постоянной скоростью перемещений, нижний пуансон жестко закреплен:
Г V | Р1 = {0,0, - V},
{ V | Р2 = {0,0,0};
2) на боковых гранях определены тривиальные статические граничные условия (также в скоростной форме) [29]:
n • о + (n • Vv • n)n • о - (Vv • n) • о = 0,
r e S3 u S4 uS5 u S6;
3) контактные условия: на торцах задано условие механического контакта, трение определяется законом Зибеля-Амонтона-Кулона. Скорости перемещения точек пуансона обозначим vp, поверхности монокристалла — vm. При этом участки поверхности геометрического контакта могут находиться в трех состояниях: зона прилипания, в которой скорости контактирующих точек поверхности равны
vp(r, t) = vm(r, t), r e S*a, зона отлипания — область поверхности, ранее находящаяся в контакте, к текущему моменту времени переходящая в свободную поверхность, определяется по условиям ann = 0,апп > 0:
о n (r, t ) = 0, r e Sf, зона скольжения:
vn(r, t ) = vm(r, t ),
v p(r, t ) - v m(r, t )
0T (r,t ) = -T f
T f =
\f a n
vP(r, t) - vm(r, t)|
при nn <Ty,
n*
re S *
ПРИ nn >Ty
где On = n • о, ann =n • о • n, от = о n
вектор
напряжений на площадке с нормалью n, его нормальная
..m __m
и касательная составляющие; vn , vT
- нормальная и
касательная к поверхности составляющая скорости точек монокристалла; Vр, Vр — нормальная и касательная к поверхности составляющая скорости точек пуансона; Т у — скорость изменения предела текучести материала в точке контакта; ц — коэффициент трения;
4) начальные условия:
(21)
5) перемещения и и напряжения о для материальной частицы определяются интегрированием (с учетом известного лагранжева закона движения г = г (^,')):
u(Ç ', t ) = } v (r (£, t), T)dT, £ '} e B (0),
0
, t) = оо (^ ) + }, T), T)dT, {£'} e B(0),
(22)
где В — область, занятая материалом монокристалла.
Для решения поставленной задачи была использована прямая модель первого типа [42] с использованием метода конечных элементов. При этом описываемый объект представляет собой отдельный монокристалл (в исходном состоянии), который аппроксимируется множеством конечных элементов. В данной задаче использованы четырехузловые симплекс-элементы с одной
точкой интегрирования; в таком случае каждый конечный элемент представляет собой кристаллит. В свою очередь, отклик данного кристаллита определяется с использованием физической теории упруговязкоплас-тичности (1)-(8). При этом скорости перемещений каждой узловой точки вычисляются в ходе решения краевой задачи, по которым вычисляются градиенты скоростей перемещений, которые являются кинематическими входными данными представленной модели кристаллита.
Четырехузловой симплекс-элемент подразумевает аппроксимацию с помощью четырех линейных функций формы:
ф(к)(х1, X2, X3) = а0к) + а(к)X1 + а(2к)х2 + «3к)х3, (23)
любое поле значений некоторой величины А интерполируется внутри конечного элемента следующим образом:
А (х1, х2, х3, t) = А(к} ^ )ф(к} (х1, х2, х3), (24)
где А(к-1 — значение величины А в узле к.
Общее количество симплекс-элементов, на которое была разбита расчетная область (объем монокристалла), составило около 164 тыс.; общее количество степеней свободы решаемой задачи — около 32 тыс. Пуансоны моделировались плоскими абсолютно твердыми пластинами, поверхность каждого аппроксимирована 150 треугольными конечными элементами.
Для решения поставленной краевой задачи методом конечных элементов использовался программный продукт собственной разработки, моделирование расчетной области, построение сетки и определение граничных условий производилось в продукте Ansys18, визуализация результатов производилась в свободно распространяемом продукте ParaView 5.2.0.
6. Оценка фрагментации объема монокристалла
Основной целью настоящей работы было выявление процессов фрагментации исходного однородного монокристалла на объемы, различающиеся эволюцией неупругой деформации и ориентациями кристаллической решетки.
На рис. 2 показано распределение интенсивности скорости пластической деформации в областях монокристалла, расположенных в окрестности грани CDD1C1, при величине относительной осадки 0.20.35 %. Из представленных результатов следует, что пластические деформации возникают в первую очередь в объемах, прилегающих к ребрам граней образца, контактирующих с пуансонами, после чего распространяются в объем монокристалла. Более подробное описание результатов исследования процесса возникновения пластических сдвигов представлено в работе [43], где основное внимание было уделено процессам, протекающим при осадке до 1 %. В [43] расчеты проводились с помощью коммерческого пакета Abaqus при помощи пользовательской процедуры отклика материала иМАТ. При этом в качестве недостатков применения иМАТ можно отметить, что отсутствовала возможность корректного описания геометрической нелинейности, а для расчета требовалось на порядок больше времени, чем при расчете с использованием программного продукта собственной разработки. В результате анизотропии и неоднородности деформирования монокристалла возникает фрагментация на 7 объемов (трех типов) с различной интенсивностью пластического деформирования, что также было предсказано теоретически и показано в натурном эксперименте [24]. В работе [43] показано, что в данных объемах активными являются различные наборы систем скольжения. Следует отме-
I
2.4520
1.8393
1.2262
0.6131
0.0000
Рис. 2. Распределение интенсивности скорости пластических деформаций в объеме монокристалла, примыкающем к грани СЬЬС1, при относительной осадке 0.20-0.35 %
I
3.112
2.334
1.556
0.778
0.000
Рис. 3. Распределение интенсивности скорости пластических деформаций в объеме монокристалла, примыкающем к грани СОВ1С1, при относительной осадке 1.5-6.0 %
тить, что в областях, примыкающих к границам указанных объемов, наблюдаются наибольшие градиенты интенсивности пластической деформации и наибольшие различия поворотов решетки.
В процессе осадки монокристалла распределение и значение интенсивности скорости пластических деформаций претерпевает несущественные изменения, что иллюстрирует рис. 3. Стоит отметить, что форма фрагментов, полученная при помощи Abaqus + иМАТ, имеет отличия, например, центральный фрагмент имел форму, близкую к правильному шестиграннику. Объем исходного монокристалла претерпевает фрагментацию на объемы с различной интенсивностью скорости неупругих деформаций, которая сохраняется на протяжении всего процесса (до 6 % относительной осадки). Вслед-
ствие наличия интенсивных сдвиговых деформаций в зоне III, материал перемещается преимущественно вдоль направления [001] [43], при этом движение материала в приторцевых объемах I затруднено вследствие возникновения существенных сил трения. Указанный характер перемещения материала приводит к значительным разворотам в областях II.
Из результатов, представленных на рис. 4, видно, что максимальное абсолютное значение угла поворота кристаллической решетки вокруг оси [110] (перпендикулярной плоскости рисунка) достигается в объемах II, а минимальное — в I и III. Ранее экспериментальным путем было установлено [44], что в областях II (на поверхности грани СВБ1С1) следы сдвига по плоскостям (111), (111) составляют угол порядка 3° с теми
ц 4.3370 3.2531
2.1687
1.0844
0.0000
Рис. 4. Распределение на грани СВБ1С1 абсолютного значения угла (в градусах) поворота кристаллической решетки в точках материала вокруг оси Х1 лабораторной системы координат при относительной осадке 1.5-6.0 %, стрелками показано направление вращения кристалличес кой решетки
же плоскостями в недеформированном состоянии. Наличие неоднородного поворота кристаллической решетки приводит к ее искривлению.
На рис. 5 представлено распределение нормы тензора кривизны л/к: кт (15) для степени относительной осадки 1.5-6.0 %. Данные результаты свидетельствуют о том, что в исходно однородном монокристалле возникают разориентированные фрагменты с достаточно четкой границей. Толщина границы не превышает 0.3 мм (высота исходного образца 6 мм). В процессе деформирования взаимное расположение фрагментов остается неизменным вплоть до 6 % относительной осадки.
На рис. 6 представлено распределение компоненты к21 тензора кривизны к (для степени относительной осадки 1.5-6.0 %), которая характеризует кривизну в плоскости, ортогональной вектору к1 подвижной системы координат, материальной линии 12 кристалличес-
кой решетки (огибающей касательных к вектору к2 подвижной системы координат). При этом максимальное отклонение вектора подвижной системы координат к1 от вектора лабораторной системы координат е1 не превышает 0.5° в областях К, следовательно, плоскость, ортогональная к1; практически совпадает с плоскостью ХХз лабораторной системы координат. Максимальное абсолютное значение кривизны достигается в областях К (рис. 6). В натурном эксперименте в данных областях наблюдается искривление следов (лежащих в плоскости XX) выхода краевых дислокаций, реализующих сдвиг по плоскостям (111), (111) (рис. 7), что является косвенным признаком искривления кристаллической решетки.
Следует отметить наличие тесной органической связи между неоднородностью пластических сдвигов и поворотами — кручениями решетки, о чем неоднократ-
0.18
0.09
0.00
-0.09
-0.18
Рис. 6. Распределение компоненты 21 кривизны решетки на грани CDDlCl монокристалла при относительной осадке 1.56.0 %, эллипсами отмечены области с максимальным абсолютным значением
[110]
0.18
0.09
0.00
0.09
0.18
Рис. 7. Рельеф грани СВБ1С1 монокристалла в эксперименте (а) и распределение компоненты к21 кривизны материала на той же грани при относительной осадке 6.0 % (б)
но отмечалось в работах В.Е. Панина и его коллег (см., например, [45]). Действительно, повороты решетки существенно зависят от пластических сдвигов по системам скольжения. Чем выше неоднородность пластических сдвигов, тем больше различие в поворотах решетки, а следовательно, тем выше абсолютные значения компонент и нормы тензора кривизны-кручения. Данное утверждение полностью согласуется с полученными теоретическими результатами: из сопоставления результатов, приведенных на рис. 2, 3 и 5, нетрудно видеть, что области с наибольшими кривизнами-кручениями совпадают с зонами наибольших градиентов пластических деформаций.
7. Выводы
Представлены результаты моделирования одноосной осадки монокристаллического образца. Для решения поставленной задачи использовалась прямая модель первого типа, основанная на методе конечных элементов. Отклик материала определялся физической упруго-вязкопластической моделью, которая явным образом учитывает сдвиги по реальным кристаллографическим плоскостям. В результате моделирования было показано, что исходный однородный объем монокристалла в процессе деформирования разделяется на 7 фрагментов, в которых интенсивность сдвиговой деформации отличается на несколько порядков, при этом максимальное значение достигается в центральной области образца, а минимальное — в приторцевых областях. Вследствие неоднородности пластических сдвигов возникает разориентация выделившихся фрагментов. Максимальное значение разориентации соседних фрагментов достигает 8°, поворот кристаллической решетки в них относительно недеформированного состояния достигает 4° при относительной осадке 6 %. В результате
неоднородного вращения кристаллической решетки в объеме образца (между фрагментами) возникают области с существенно искривленной решеткой. В численном эксперименте получено, что в ходе осадки в ряде областей происходит искривление кристаллической решетки монокристалла, что согласуется с результатами натурного эксперимента, в котором показано, что на поверхности также наблюдаются искривленные линии скольжения краевых дислокаций.
Работа выполнена при поддержке гранта РНФ 1719-01292.
Литература
1. Теплякова Л.А., Куницына Т.С., Конева Н.А. Влияние ориентации монокристаллов сплава со сверхструктурой L^rn закономерности формирования субструктуры при сжатии // Деформация и разрушение материалов. - 2015. - № 2. - С. 40-48.
2. Теплякова Л.А., Беспалова И.В., Лычагин Д.В. Пространственная организация деформации в [1-12]-монокристаллах алюминия при сжатии // Физ. мезомех. - 2009. - Т. 12. - № 2. - С. 67-76.
3. Lychagin D.V., Tarasov S.Yu., Chumaevskii A.V., Alfyorova E.A. Macro-
segmentation and strain hardening stages in copper single crystals under compression // Int. J. Plasticity. - 2015. - V. 69. - P. 36-53.
4. Saai A., Louche H., Tabourot L., Chang H.J. Experimental and numerical study of the thermo-mechanical behavior of Al bi-crystal in tension using full field measurements and micromechanical modeling // Mech. Mater. - 2010. - V. 42. - P. 275-292.
5. Groh S., Marin E.B., Horstemeyer M.F., Zbib H.M. Multiscale modeling of the plasticity in an aluminum single crystal // Int. J. Plasticity. -2009. - V. 25. - P. 1456-1473.
6. Dumoulin S., Tabourot L. Experimental data on aluminium single crystal behaviour // Proc. Inst. Mech. Eng. C. - 2005. - V. 219. - No. 9. -P. 1159-1167.
7. Liu M., Lu C., Tieu A.K. Crystal plasticity finite element method modelling of indentation size effect // Int. J. Solid. Struct. - 2015. - V. 54. -P. 42-49.
8. Gerard C., Cailletaud G., Bacroix B. Modeling of latent hardening produced by complex loading paths in fcc alloys // Int. J. Plasticity. -2013. - V. 42. - P. 194-212.
9. Wang Z.Q., Beyerlein I.J., LeSar R. Plastic anisotropy in fcc single crystal in high rate deformation // Int. J. Plasticity. - 2009. - V. 25. -P. 26-48.
10. Preussner J., Rudnik Y., Brehm H., VolklR., Glatzel U. A dislocation density based material model to simulate the anisotropic creep behavior of single-phase and two-phase single crystals // Int. J. Plasticity. -2008. - doi 10.1016/j.ijplas.2008.04.006.
11. Ha S., Jang J.-H., Kim K.T. Finite element implementation of dislocation-density-based crystal plasticity model and its application to pure aluminum crystalline materials // Int. J. Mech. Sci. - 2017. -V. 120. - P. 249-262. - doi 10.1016/j.ijmecsci. 2016.11.011.
12. Harder J. FEM-simulation of the hardening behavior of FCC single crystals // Acta Mech. - 2001. - V. 150 - Р. 197-217.
13. Horstemeyer M.F., Potirniche G.P., Marin E.B. Crystal Plasticity // Handbook of Materials Modeling / Ed. by S. Yip. - Netherlands: Springer, 2005. - Р. 1133-1149.
14. Roters F. Application of crystal plasticity FEM from single crystal to bulk polycrystal // Comput. Mater. Sci. - 2005. - V. 32. - Р. 509-517.
15. Staroselsky A., Anand L. A constitutive model for hcp materials deforming by slip and twinning: application to magnesium alloy AZ31B // Int. J. Plasticity. - 2003. - V. 19. - Р. 1843-1864.
16. Wang H., Wu P.D., Wang J., Tome C.N. A crystal plasticity model for hexagonal close packed (HCP) crystals including twinning and de-twinning mechanisms // Int. J. Plasticity. - 2013. - V. 49. - Р. 3652.- doi 10.1016/j .ijplas.2013.02.016.
17. Zaafarani N., Raabe D., Singh R.N., Roters F., Zaefferer S. Three-dimensional investigation of the texture and microstructure below a nanoindent in a Cu single crystal using 3D EBSD and crystal plasticity finite element simulations // Acta Mater. - 2006. - V. 54. - P. 1863— 1876.
18. Хирт Дж., Лоте И. Теория дислокаций. - М.: Атомиздат, 1972. -600 с.
19. Хоникомб Р. Пластическая деформация металлов. - М.: Мир, 1972. - 408 с.
20. Asaro R.J. Geometrical effects in the inhomogeneous deformation of ductile single crystals // Acta Metall. - 1979. - V. 27. - No. 3. - P. 445453.
21. Peirce D., Asaro R., Needleman A. An analysis of nonuniform and localized deformation in ductile single crystals // Acta Metall. - 1982. -V. 30. - P. 1087-1119.
22. Szczerba M., Korbel A. Strain softening and instability of plastic flow in Cu-Al single crystals // Acta Metall. - 1987. - V. 35. - No. 5. -P. 1129-1135.
23. Теплякова Л.А., Лычагин Д.В., Козлов Э.В. Локализация сдвига при деформации монокристаллов алюминия с ориентацией оси сжатия [001] // Физ. мезомех. - 2002. - Т. 5. - № 6. - С. 49-55.
24. Теплякова Л.А., Лычагин Д.В., Беспалова И.В. Закономерности макролокализации деформации в монокристаллах алюминия с ориентацией оси сжатия [110] // Физ. мезомех. - 2004. - Т. 7. -№ 6. - С. 63-78.
25. Chang Y.W., Asaro R.J. An experimental study of shear localization in aluminum-copper single crystals // Acta Metall. - 1981. - V. 29. -No. 1. - P. 241-257.
26. Фридель Ж. Дислокации. - М.: Мир, 1967. - 644 с.
27. Миркин Л.И. Физические основы прочности и пластичности (Введение в теорию дислокаций). - М.: Изд-во МГУ, 1968. - 538 с.
28. Новиков И.И. Дефекты кристаллического строения металлов: учебное пособие. - М.: Металлургия, 1975. - 208 с.
29. Поздеев А.А., Трусов П.В., Няшин Ю.И. Большие упругопласти-ческие деформации: теория, алгоритмы, приложения. - М.: Наука, 1986. - 232 с.
30. Трусов П.В., Швейкин А.И., Янц А.Ю. О разложении движения, независимых от выбора системы отсчета производных и определяющих соотношениях при больших градиентах перемещений: взгляд с позиций многоуровневого моделирования // Физ. мезо-мех. - 2016. - Т. 19. - № 2. - С. 47-65.
31. Панин В.Е., Егорушкнн В.Е., Панин А.В. Физическая мезомеха-ника деформируемого твердого тела как многоуровневой системы. I. Физические основы многоуровневого подхода // Физ. мезомех. -2006. - Т. 9. - № 3. - С. 9-22.
32. Макаров П.В. Математическая теория эволюции нагружаемых твердых тел и сред // Физ. мезомех. - 2008. - Т. 11. - 3. - С. 1935.
33. ТрусовП.В., Швейкин А.И. О разложении движения и определяющих соотношениях в геометрически нелинейной упруговязкоплас-тичности кристаллитов // Физ. мезомех. - 2016. - Т. 19. - № 3. -С. 25-38.
34. Трусов П.В., Швейкин А.И., Нечаева Е.С., Волегов П.С. Многоуровневые модели неупругого деформирования материалов и их применение для описания эволюции внутренней структуры // Физ. мезомех. - 2012. - Т. 15. - № 1. - С. 33-56.
35. Трусделл К. Первоначальный курс рациональной механики сплошных сред. - М.: Мир, 1975. - 592 с.
36. Трусов П.В., Нечаева Е.С., Швейкин А.И. Применение несимметричных мер напряженного и деформированного состояния при построении многоуровневых конститутивных моделей материалов // Физ. мезомех. - 2013. - Т. 16. - № 2. - С. 15-31.
37. Трусов П.В., Янц А.Ю. О физическом смысле неголономной меры деформации // Физ. мезомех. - 2015. - Т. 18. - № 2. - С. 13-21.
38. Asaro R.J., Needleman A. Texture development and strain hardening in rate dependent polycrystals // Acta Metall. - 1985. - V. 33. - No. 6. -P. 923-953.
39. Янц А.Ю. Двухуровневая модель для описания неупругого деформирования поликристаллов: приложение к анализу сложного на-гружения в случае больших градиентов перемещений: Дис. ... канд. физ.-мат. наук. - Пермь: ИМСС УрО РАН, 2016.
40. Pietraszkiewicz W., Eremeyev V.A. On natural strain measures of the nonlinear micropolar continuum // Int. J. Solid. Struct. - 2009. -V. 46. - No. 3-4. - P. 774-787.
41. Гордеев В.Н. Кватернионы и бикватернионы с приложениями в геометрии и механике. - Киев: Сталь, 2016. - 316 с.
42. Трусов П.В., Швейкин А.И. Многоуровневые физические модели моно- и поликристаллов. Прямые модели // Физ. мезомех. - 2011. -Т. 14. - № 5. - С. 5-30.
43. Янц А.Ю., Теплякова Л.А. Моделирование нагружения монокристаллических образцов // Изв. вузов. Физика. - 2015. - Т. 58. -№ 3. - С. 330-335.
44. Беспалова И.В. Влияние кристаллогеометрических характеристик на закономерности фрагментации и локализации сдвиговой деформации в монокристаллах алюминия при сжатии // Дис. ... канд. физ.-мат. наук. - Томск: ИФПМ СО РАН, 2008.
45. Panin V.E., Egorushkin V.E., Surikova N.S., Pochivalov Yu.I. Shear bands as translation-rotation mode of plastic deformation in solids under alternate bending // Mater. Sci. Eng. A. - 2017. - V. 703. -P. 451-460.
Поступила в редакцию 14.11.2017 г., после переработки 28.12.2017 г.
Сведения об авторах
Трусов Петр Валентинович, д.ф.-м.н., проф., зав. каф. ПНИПУ, tpv@matmod.pstu.ac.ru Янц Антон Юрьевич, ассист. ПНИПУ, maximus5.59@gmail.com Теплякова Людмила Алексеевна, д.ф.-м.н., проф. ТГАСУ, lat168@mail.ru
