Несимметричная физическая теория пластичности ГЦК-поликристаллов: особенности численной реализации некоторых схем деформирования

УДК 539.3

П.С. Волегов, А.Ю. Янц

Пермский государственный технический университет

НЕСИММЕТРИЧНАЯ ФИЗИЧЕСКАЯ ТЕОРИЯ ПЛАСТИЧНОСТИ ГЦК-ПОЛИКРИСТАЛЛОВ: ОСОБЕННОСТИ ЧИСЛЕННОЙ РЕАЛИЗАЦИИ НЕКОТОРЫХ СХЕМ ДЕФОРМИРОВАНИЯ*

В статье рассмотрены вопросы, связанные с построением физической теории пластичности поликристалла, учитывающей развороты кристаллических решеток зерен, на базе несимметричных мер напряженного и деформированного состояния. Предложена двухуровневая модель упруговязкопластического деформирования поликристаллического агрегата, записан общий вид несимметричного закона упругости в скоростной релаксационной форме. Сформулированы соотношения для поверхностных и объемных моментов, приводящих к разворотам кристаллической решетки зерен и фрагментов. Отдельно рассмотрены проблемы, возникающие при численной реализации некоторых часто используемых схем деформирования - осадки и стесненной осадки. Показано, что выбор неявной схемы интегрирования второго порядка точности позволяет существенно ускорить процесс численного счета за счет возможности увеличения шага интегрирования при сохранении требуемой точности.

Ключевые слова: физические теории, несимметричная упруговязкопластичность, ротация, текстура, численные методы, схема Адамса-Мултона.

Введение

В физических теориях пластичности, построенных на основе симметричных мер напряженного и деформированного состояния, существуют проблемы, связанные с неоднозначностью определения активных систем скольжения дислокаций, а также с невыполнением уравнения баланса момента количества движения при рассмотрении процессов, связанных с глубокой перестройкой мезоструктуры поликристаллов; при этом отказываться от учета распределенных поверхностных моментов, приводящих к появлению моментных напряжений, в задачах эволюции микро- и мезоструктуры нельзя, так как это приведет к невозможности физически прозрачно описывать, например, повороты кристаллитов, а следовательно, текстурирование.

*Работа выполнена при финансовой поддержке РФФИ (гранты 10-08-96010-р_урал_а, 10-08-00156-а).

Целью работы является построение математической модели неупругого деформирования ГЦК-поликристаллов на базе физической теории пластичности, в которой не используется симметризация мер напряженного и деформированного состояний, а также их скоростей.

В работах [1-3] предлагается подход к построению многоуровневых моделей (в частном случае - двухуровневых) деформируемого по-ликристаллического агрегата (ПКА), в частности общая структура конститутивных моделей с внутренними переменными для описания процессов эволюции микро- и мезоструктуры материала. В рамках данного исследования также будет использоваться двухуровневая модель упругопластического деформирования ПКА. Для простоты понимания величины, относящиеся к макроуровню, обозначаются заглавными буквами, родственные параметры, определенные на мезоуровне - теми же строчными. Задается схема деформирования на макроуровне, для деформированного состояния на мезоуровне (зерно, субзерно) принимается модифицированная гипотеза Фойгта для меры скорости деформации. В роли определяющего соотношения на макроуровне используется закон Г ука в скоростной релаксационной форме; скорость пластических деформаций определяется из модели второго масштабного уровня (мезоуровня) по скоростям сдвигов по активным системам скольжения, для поиска которых, в свою очередь, используется определяющее соотношение вязкого типа. Предполагается, что процессы деформирования являются квази-статическими и протекают при низких гомологических температурах, что позволяет не учитывать диффузные механизмы (рекристаллизацию, возврат, образование атмосфер примесных атомов (Коттрелла, Сузуки).

1. Несимметричные меры скорости деформации и напряжения

Для описания процессов деформирования сплошной среды предлагается использовать в качестве меры скорости деформации транспонированный градиент скорости [4]:

? - V/.

Для градиента места используется мультипликативное разложение на обратимую (упругую) и необратимую части: ї = їе  їт . Подставив это соотношение в выражение для £ , получаем

£ __ е ї іп (ї іп )-1 (ї е )-1 ^ — е ї іп + ї е — іп ^ (ї іп )-1 ^ ї е ^ 1

= і-е  їіп

+ їе ї ї )-1 |(їе )-1.

Первое слагаемое суммы описывает обратимые деформации среды, второе, вообще говоря, в записанном виде не может быть названо скоростью необратимых деформаций, для этого необходимо перейти к разгруженной конфигурации.

Далее введем на мезоуровне меру деформированного состояния q, формально используя определение меры скорости деформации как производной от соответствующей (неголономной) меры деформации:

qr _ <1 - ю  q + q  ю _ £ _ уУ .

В качестве спина на мезоуровне предлагается использовать спин решетки ю.

2. Несимметричный закон упругости мезоуровня

В классической линейной механике закон упругости получается из разложения в ряд свободной энергии в окрестности нулевых упругих деформаций; с учетом начального ненапряженного и недеформи-рованного состояния в итоге получается соотношение вида

д 2Б

°у _ р д е д е чеи _с укіїм; _с: іє ,

дЧ у дЧ к1

где с - тензор четвертого ранга, описывающий упругие свойства материала, р - плотность, Б - свободная энергия Гельмгольца. Из вида (4) можно утверждать симметрию по парам индексов, сі]к1 _ сИі]. Утверждать симметрию с внутри пар индексов нельзя в силу несимметрии на мезоуровне мер напряженного и деформированного состояний.

В [5] показано, что вместе с условием симметрии по парам индексов такой алгоритм позволяет утверждать, что для материала с кубической симметрией решетки тензор упругих свойств с всего имеет 4 независимых компоненты - с1111, с1122, с1212, с1221 (в кристаллографической системе координат).

Тогда на основании вышесказанного несимметричный закон упругости на мезоуровне можно записать в виде

о _ с: ,

здесь с - тензор четвертого ранга, имеющий 4 независимые компоненты,

или в скоростной форме

ог _ с: £е,

г

где о - соответствующая коротационная производная мезоуровня, введенная выше и учитывающая вращение кристаллографической системы координат, связанной с решеткой, как жесткого целого. Приведем также закон упругости макроуровня:

Яг = С:

где г и Zп - меры скорости деформации и скорости необратимых

деформаций на макроуровне соответственно.

Для решения задачи идентификации вновь появившихся компонент тензора упругих свойств потребовалось поставить и реализовать серию численных экспериментов на простое нагружение, позволивших определить значения всех независимых компонент. В статье [5] при помощи моделирования методом молекулярной статики проводится качественный анализ поведения диаграммы напряжение-деформация для монокристалла при различных отклонениях компонент с1212 , с1221

тензора упругих характеристик друг от друга. Показано, что даже незначительное расхождение этих компонент тензора упругих свойств существенным образом влияет на поведение представительного объема материала.

3. Описание ротации решетки и фрагментации зерен

Основной причиной разворотов решетки зерен будем считать несовместность сдвигов по системам скольжения в соседних зернах (моделирующих, в свою очередь, движение дислокаций). Более подробно проблема описания ротаций кристаллической решетки рассмотрена в работе [4], здесь остановимся на этом вопросе кратко. Скорость изменения вектора поверхностного момента, действующего на часть границы анализируемого зерна (фрагмента зерна) в результате сопротив-

ления переходу дислокаций из анализируемого зерна (фрагмента) в соседние (т = 1,...,М ), можно определить как сумму:

л

Ъ

п=\

X N х

'V

где (-)г - соответствующая коротационная производная [4], X - экспериментально определяемый (в Па • м ) параметр, N - внешняя для анализируемого фрагмента единичная нормаль к границе с соседним фрагментом, у ', у](т) - скорости сдвигов, Ь!, Ь](т) - единичные векторы по направлениям векторов Бюргерса, п', п](т) - нормали для систем скольжения в исследуемом и соседнем фрагментах соответственно, К - число систем скольжения (для ГЦК-кристалла с учетом удвоения - 24).

Для того чтобы определять величину спина решетки элемента ротации (ЭР), требуется «свести» все динамические причины разворота в некоторую эквивалентную величину, в качестве которой предлагается использовать удельный объемный момент:

\ М

м' = Т7ТЪ[тХ + гт х

* т=1

где М ' ([м! ^ = Па) - вектор объемного момента, действующий на данный ЭР ', Бгт - площадь плоского участка границы т , разделяющей два

ЭР, V' - объем данного ЭР, гт - радиус-вектор, проведенный от центра масс ЭР к средней точке фасетки т , & - тензор напряжений Коши, действующих в данной подобласти ЭР, Nт - внешняя нормаль к фасетке т . Второе слагаемое правой части (9) определяет дополнительную составляющую момента, возникающую на границе ЭР вследствие несимметричного (в общем случае) напряженного состояния.

Определив элементы ротации, необходимо рассчитать скорость ротации кристаллографической системы координат данного ЭР. Описание разворота на каждом шаге интегрирования производится при помощи ортогонального тензора ДR = (со8Дф +1^ + С08Дф! +

+8тДф e х I (I - единичный тензор второго ранга), определяющего поворот системы координат вокруг мгновенной оси вращения e на некоторый угол Дф. При этом

• направление e считается соосным вектору объемного момента:

M

e = -

• для величины угловой скорости разворота принимается гипотеза об аддитивности и гипотеза единой кривой:

(р -

— №1 + — IM|, |м| - Мг и M • М > 0,

А н сю)

А |М|, в противном случае,

где А, Н - экспериментально определяемые параметры материала.

Критическое значение Мс в первом приближении будем вычислять следующим образом:

1 № п 1

Мс = 17Ег' х

V %-I1-

где Мс, Па - критическое значение объемного момента, Sg - площадь фасетки входящего в рассматриваемый ЭР граничного фрагмента g, тg, Па - среднее критическое напряжение по системам скольжения граничного элемента, Ng - количество граничных фрагментов данного

ЭР, V - объем элемента ротации; суммирование проводится по «граничным» фрагментам ЭР.

4. Проблема численного описания одноосного напряженного состояния на макроуровне

В натурном эксперименте при одноосном напряженном состоянии имеется одна ось, вдоль которой производится нагружение: при жестком - задаются перемещения вдоль данной оси, при мягком - усилия. Используемая в данной работе двухуровневая модель основана на гипотезе Фойгта, служит для связи переменных мезо- и макроуровней, что требует задания только жесткого нагружения. По условиям одноосного нагружения существуют две свободные поверхности, перемещения которых в общем случае неизвестны, отсюда вытекают некото-

рые особенности численного моделирования для реализации одноосного напряженного состояния на макроуровне.

Положим, что деформирование производится вдоль оси ОХ3 лабораторной системы координат, тогда при жестком нагружении задается фиксированное перемещение торцов образца вдоль оси ОХ3. В силу того что задача решается в скоростях, имеем:

Z33 = а - фиксированное;

Я(>. = 0, і, ] Ф 3, Я33 Ф 0 - условие одноосности;

С - тензор упругих характеристик на макроуровне;

Яг = С: ^-Zш), Яг = Я-П• Я + Я• П, закон Гука на макроуровне;

Zш - осредненные пластические мезодеформации.

Из системы (12) необходимо получить оставшиеся 5 неизвестных компонент Z и компоненту Я33. Используем закон Гука - неизвестные оставим слева, известные величины перенесем вправо:

'Я = С: ^ - Ziи) + П • Я - Я • П,

Я - С^ = -С^іи + П • Я - Я • П.

Отсюда, расписывая покомпонентно, учитывая (12) и факт симметричности всех характеристик на макроуровне, получаем систему линейных алгебраических уравнений (СЛАУ) в обычном матричном виде:

-1 с С3311 2С 2С3312 2с 2С3313 с с3322 2С ^ 2С3323 (я ^ ¿-,33

0 с с1111 2С 2С1112 2с 2С1113 с с1122 2С 2С1123

0 с с1211 2С 2С1212 2с 2С1213 с с1222 2С 2С1223 712

0 с с1311 2с 2С1312 2с 2С1313 с с1322 2С 2С1323 7 713

0 с с2211 2с 2С2212 2с 2С2213 с с2222 2С 2С2223 7 7 22

0 с с2311 2с 2С2312 2с 2С2313 с с2322 2С 2С2323 У 7 V 723 у

= Ь

где

ь =

С3333^33 + I I С33у^;, ¿=1 ; =1

-С, ,33 г33 + 1 I ;

¿=1 _/=1

33

-с^+¿1 е,гяг;

1=1 ;=1

33

-С^ +11С01Д + («31 - П.3>2!!

¿=1 7=1

33

-С2!33233 +11 ад”

¿=1 7=1

-ад„+¿1 ад;+(о,, - п!3>Е33

¿=1 ; =1

Данная СЛАУ разрешается на каждом шаге интегрирования для нахождения скоростей деформаций и компоненты Ё33, согласующихся

с условием одноосности.

Ниже (рис. 1, 2) приведены зависимости компонент тензора напряжений £ на макроуровне от интенсивности накопленных деформаций при использовании численного описания одноосного напряженного состояния, представленного выше, и при использовании псевдоод-ноосного НДС с тензором скоростей деформаций (13), который часто используется в работах для описания одноосного напряженного состояния:

Ъи

а/ /2 0 0

0 а 2 0

0 0 -а

Рис. 1. Компоненты напряжений для «псевдоодноосного НДС»

Рис. 2. Компоненты напряжений для одноосного НДС

Из графиков видно, что в сравниваемых двух случаях имеем существенно различные напряженно-деформированные состояния, предел текучести также отличается на ~20% (одноосное НДС - 41 МПа, «псевдоодноосное» - 33 МПа), что свидетельствует о невозможности идентификации модели на натурных одноосных экспериментах с применением в численных расчетах деформации вида Х*итса .

5. Численная реализация стесненной осадки

Эксперимент на стесненную осадку производится аналогично эксперименту по одноосному сжатию, за исключением того, что перемещения свободных поверхностей вдоль одной из осей (ОХ\) невозможны, причем трение на этих поверхностях отсутствует. В этом случае получаем четыре неизвестных компоненты Z и две компоненты скорости напряжений. Имеем для сжатия вдоль ОХ3 при стеснении вдоль ОХі:

Zзз = а - фиксированное;

Аі = 0;

Zш - осредненные пластические мезодеформации;

Еі/ = 0, Е33 ^ 0,Е 11 ^ 0;

С - тензор упругих характеристик на макроуровне;

Ег = С: ^ - Zш ), Ег = Е - П • Е + Е • П, закон Гука на макроуровне.

Аналогично схеме, приведенной для одноосного сжатия, имеем:

Е = С: ^ - Z,'и) + П • Е - Е • П,

Е - С : Z = -С^'и + П • Е - Е • П.

Отсюда, расписывая покомпонентно, учитывая (14) и факт симметричности всех характеристик на макроуровне, получаем систему линейных алгебраических уравнений (СЛАУ) в обычном матричном виде:

-1 0 2С 2С3312 2С 2С3313 С С3322 2С ^ 2С3323 (Ё ^ 33

0 -1 2С 2С1112 2С 2С1113 С С1122 2С 2С1123 Ё 11

0 0 2С 2С1212 2С 2С1213 С С1222 2С 2С1223 ¿12

0 0 2С 2С1312 2С 2С1313 С С1322 2С 2С1323 7 ¿13

0 0 2С 2С2212 2С 2С2213 С С2222 2С 2С2223 7 7 22

0 0 2С 2С2312 2С 2С2313 С С2322 2С 2С2323 У 7 V723 у

= Ь

где

ь =

-С3333 2зз+11С33,2,

,=1 ,=1

-С„зз2зз+1 1 Сщг-;

,=1 , =1

-с^А +11 С1г ,2- + (п„ - П21)ЕП

,=1 , =1

-СШз2зз + 11 Сщг; +

,=1 , =1

33

-С2233233 +11С22,21

?=1 У =1

33

С2ззз2зз + 11 СШ12‘; + (П32 -Й2з)Е

?=1 У =1

33

Данная СЛАУ разрешается на каждом шаге интегрирования для нахождения необходимых для выполнения условий (14), скоростей деформаций и компонент 133,111.

Ниже приведены зависимости компонент тензора напряжений Е на макроуровне от интенсивности накопленных деформаций при использовании численного описания стесненной осадки, представленной выше, и при использовании псевдоодноосного НДС с тензором скоростей деформаций.

Компоненты напряжений, МПа 0 -5-1 0-15-20-25-30 -40- -50- V , '•*

1 ,^

--4-4-4--

—

! ! I

-60

1 33 ...

Интенсивность накопленных деформаций, %

Рис. 3. Компоненты напряжений для стесненной осадки, определенные согласно

Из рис. 3 видно, что нетривиальными являются лишь компонента Е33 - напряжения вдоль оси ОХ3 и компонента Е11 - напряжения вдоль

оси стеснения, что свидетельствует о приемлемости вышеизложенных соотношений и рассуждений.

Зачастую в работах при описании стесненной осадки используется тензор скоростей деформации вида

Ъс

0 0 0

0 а 0

0 0 -а

Из рис. 3-4 нетрудно заметить, что компоненты макронапряжений при строгом задании стесненной осадки на макроуровне и при задании деформирования в виде Ъсот1г значительно отличаются, что говорит о совершенно разных НДС в этих случаях.

-

1 22

2

11

у1

213 2 12 ^23

2

33

О 2 4 6 8 10 12 14 16 18 20 22 24 26 28

Интенсивность накопленных деформаций, %

Рис. 4. Компоненты напряжений для стесненной осадки в приближении с тензором

1 гу сотії

скорости деформации на макроуровне вида Л

При этом предел текучести отличается на ~30 % (строгое задание стесненной осадки НДС - 52МПа, при использовании х,соШг -35 МПа), что свидетельствует о невозможности идентификации модели на натурных экспериментах при стесненной осадке с применением в численных расчетах деформации вида %соП!1‘г.

6. Схема численного интегрирования Адамса-Мултона

Схема Адамса-Мултона является неявной схемой интегрирования обыкновенных дифференциальных уравнений (ОДУ). Известно, что неявные схемы обладают большей устойчивостью по сравнению с явными схемами; при этом схема Адамса-Мултона хорошо подходит для ОДУ с сильно нелинейными правыми частями, что имеет место в рассматриваемом случае.

Пусть имеем ОДУ вида у = /(х, у) , тогда схема интегрирования Адамса-Мултона имеет следующий вид:

таким образом, сразу виден недостаток данной схемы - неявное разрешение относительно значения переменной на следующем шаге интегрирования. Разрешение такого нелинейного уравнения сводится к двухшаговой итерационной процедуре:

Схема Адамса-Мултона была использована при решении рассматриваемых задач, что позволило повысить точность вычислений, а также значительно увеличить скорость счета по сравнению с явной схемой Эйлера, несмотря на дополнительные итерации для решения системы нелинейных уравнений и хранение дополнительной информации - значений правых частей с предыдущих итераций. Ниже для наглядности приведен общий алгоритм (без подробного рассмотрения упрочнения и ротации) применительно к задаче определения НДС поликристалла при произвольном жестком нагружении.

Решение задачи на каждой итерации проводится полностью (где это возможно) в скоростях и потом интегрируется.

1. Пусть имеем некоторый фрагмент (номер фрагмента для упрощения записи ниже опустим) на шаге к, для которого известны напряжения на начало шага а(к), скорость напряжений а(к), сдвиги по СС у(к), і —1,24

скорости сдвигов по СС у(к), і — 1,24, критические напряжения

т.к), і — 1,24, расчет для каждого фрагмента проводится в КСК.

2. С помощью явной схемы Эйлера выполняем первую итерацию: а) находим критические напряжения в СС:

г

Мі)

— ) +

/ (у (‘У(к)) Д», і —1,24;

б) находим скорости сдвигов и сдвиги в СС [6]:

У

у 0 Н ('

.

г

24

1/т

“8”(г(к)),т(к) — Ь,.п,: ч,к),

¿,(к) =^у (к тЬп(,

і—1

у( к+1)(1) — у(к) +у(к )(1) Д», і —1,24; в) определяем напряжения:

,

' — С :(Л - Ліп(к)),

іи

(к )(1) Д».

+ <т

3. Выполняем вторую итерацию схемы Адамса-Мултона: а) определяем критические напряжения:

г(к+1)(2) — г(к) + Іс(і) Іс(і) ^

Д»

Т

/ (у'*). У (‘)) + / (у (к+И1). У (‘)(1))]. і —1.24;

б) находим скорости сдвигов в СС:

У і

:у0Н (г(

с

24

1/т

Vм) —Еу;‘ )(!)ЬіПі,

і—1

Уі

- У(‘1+Д\ [У “ ^ +У (к >(2> ],

У(к )(0) —У(к-1), і —1,24.

4. Выполняем итерацию для уточнения напряжений, значение напряжений на конец данного шага и на начало следующего - о(к+1)(2):

О(к )(2) = с: (Ъ - Ъп(к)),

0

2 1 О

5. На конец шага пересчитываем ориентацию фрагмента по явной схеме Эйлера, так как ориентации кристаллитов влияют только на макрохарактеристики поликристаллического агрегата и полностью определяются мезонапряжениями и неупругими мезодеформациями:

м(к) = - У

vm¿

\Х № X

24

Уг

п,Ь,

24

У (к

+ гт X

(Nт • о(к+1)(2)) БШ

м

ф(к) =

м

—М

А

—М

.А Н

< м.

= м.

1?.(к) = (С08ф(к) + 1)ее + С08ф(к )1 + 8тф(к) е XI, ю - Я(к) • Яг,

Я(к+1) - Я(к) • Я(к)А?.

6. Определяются макронапряжения:

(8уш (о(к+1)(2) )}.

Е(к+1) - (8

На рис. 5 представлена зависимость отклонения макронапряжений от некоторого решения, которое считается «точным» (в силу выбора очень маленького шага) при использовании той или иной схемы интегрирования. Видно, что схема численного интегрирования Адамса-Мултона имеет меньшую погрешность по сравнению со схемой Эйлера. Также можно отметить, что схема Эйлера пригодна для расчета, т.е. устойчива при шаге интегрирования менее 0,1 с, в то же время схема Адамса-Мултона устойчива вплоть до шага интегрирования 0,17 с. В расчетах приемлемой точностью можно считать ~1 %, поэтому шаг в данной работе взят 0,1 с, при этом скорость счета в 1,5 раза выше, чем при использовании схемы Эйлера при той же погрешности, т.е. при шаге 0,05 с. Отсюда следует целесообразность использования пред-

ставленной неявной схемы численного интегрирования Адамса-Мултона, что дает прирост в скорости нахождения «точных» расчетов, а также возможность быстрого получения «неточных» (погрешность порядка ~2 %) на больших шагах интегрирования.

2.5 2

1.5 1

0,5 0

0,02 0,03 0,04 0,05 0,07 0,08 0,09 0,1 0,12 0,14 0,16 0,17

Рис. 5. Зависимость погрешности от шага интегрирования для различных схем.

Ось абсцисс - шаг по времени

Заключение

Таким образом, в статье рассмотрены вопросы, связанные с построением физической теории пластичности поликристалла, учитывающей развороты кристаллических решеток зерен, на базе несимметричных мер напряженного и деформированного состояния. Предложена двухуровневая модель упруговязкопластического деформирования поликристаллического агрегата, записан общий вид несимметричного закона упругости в скоростной релаксационной форме. Сформулированы соотношения для поверхностных и объемных моментов, приводящих к разворотам кристаллической решетки зерен и фрагментов. Отдельно рассмотрены проблемы, возникающие при численной реализации некоторых часто используемых схем деформирования - осадки и стесненной осадки. Выбор неявной схемы интегрирования второго порядка точности позволяет существенно ускорить процесс численного счета за счет возможности увеличения шага интегрирования при сохранении требуемой точности.

Библиографический список

1. Constitutive relations and their application to the description of microstructure evolution / P.V. Trusov, V.N. Ashikhmin, P.S. Volegov, A.I. Shveykin // Physical Mesomechanics. - 2010. - Vol. 13, Is. 1-2, Janu-ary-April. - P. 38-46.

2. Трусов П.В., Волегов П.С. Определяющие соотношения с внутренними переменными и их применение для описания упрочнения в монокристаллах // Физическая мезомеханика. - Томск, 2009. - Т. 12, № 5. - С. 65-72.

3. Трусов П.В., Ашихмин В.Н., Волегов П.С. Двумерная модель пластического деформирования монокристалла // Математические модели и методы механики сплошных сред: сб. науч. ст. / Ин-т автоматики и процессов управления ДВО РАН. - Владивосток, 2007. - С. 259-269.

4. Трусов П.В., Волегов П.С., Янц А.Ю. Несимметричная физическая теория пластичности для описания эволюции микроструктуры поликристаллов // Физическая мезомеханика. - Томск, 2011. - Т. 14, № 1 -С. 19-31.

5. Волегов П.С., Шулепов А.В. Упругие константы монокристалла в несимметричной физической теории пластичности // Вестник ПГТУ. Механика. - Пермь: Изд-во Перм. гос. техн. ун-та, 2010. - № 1. -С.19-34.

6. Anand L. Single-crystal elasto-viscoplasticity: application to texture evolution in polycrystalline metals at large strains // Comput. Methods Appl. Mech. Engrg. -2004. - Vol.193. -P. 5359-5383.

Получено 21.03.2011