ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ ДЕФОРМАЦИИ МАТЕРИАЛОВ С УЧЕТОМ НЕУСТОЙЧИВОЙ ВЕТВИ σ-ε-ДИАГРАММЫ Текст научной статьи по специальности «Физика»

Численное моделирование деформации материалов с учетом неустойчивой ветви а-е-диаграммы

О.И. Черепанов

Институт физики прочности и материаловедения СО РАН, Томск, 634021, Россия

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

1. Введение

Идеально упругопластическая модель среды, в рамках которой получен ряд важных результатов теории пластичности, тем не менее, не учитывает многие эффекты в поведении реальных материалов. Для описания деформационного упрочнения материалов, связанного с многообразием дефектных структур на микро- и мезо-уровнях, моделирования эффектов локализации деформаций на неоднородностях структуры, а также изучения пластической деформации как процесса потери материалом сдвиговой устойчивости в рамках развиваемого академиком Паниным В.Е. подхода физической мезо-механики [1] необходимо применение более сложных моделей среды. Глобальная потеря сдвиговой устойчивости приводит либо непосредственно к разрушению материала, либо к предварительному появлению нисходящей ветви на диаграмме нагружения. Разупрочнение может сопровождаться аномальным увеличением объема материалов при сжатии. В опытах Бриджмена П. [2] увеличение объема при сжатии наблюдалось для мрамора. ЛатынинаЛ.А. [3] приводит данные о подобных результатах для стали, мыльного камня, талька, диабаза, гранита, объясняя эффект образованием микротрещин. Аналогичные данные для поликристаллических металлов приводит Кафка В. [4]. Математическая модель пластичности, способная описать такие эффекты, должна учитывать влияние гидростатического напряжения. Далее рассматривается расчетная схема, построенная на основе модели упругопластической среды этого класса, предложенной в работе ДрагонаА., МрузаЗ. [5]. Модель описывает эффекты упрочнения и накопления микро- или мезоповреждений материала, разупрочнение материала на этапе неустойчивого деформирования, а

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

2. Описание модели среды

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

1 т

Фц = 777-X ’' ЬГ ■ й1), (1

V к=1

где п( к ’ — вектор единичной нормали к элементарной площадке к) на поверхности трещины; ’ — разрыв смещений на ее поверхности, который характеризует взаимное проскальзывание и раскрытие берегов трещины; Уг — характерный представительный объем материала.

© Черепанов О.И., 1999

Таким образом, рассматривается континуум с плотным распределением микротрещин в локальных объемах материала Гг. В работе [5] отмечается, что при произвольной ориентации трещин в объеме Vг определить параметр разрушения фу количественно трудно. Но задача облегчается для случаев, когда обоснованы представления о существовании преимущественной ориентации (трещины параллельны) или для однородного распределения трещин по всем направлениям. В результате принятия этих допущений может быть построена модель, по форме близкая к концепции скольжения в теории пластичности [7], или же дилатанси-онная модель, учитывающая пластическое изменение объема («разрыхление» материала), накопление микро-или мезоповреждений среды в процессе деформирования. В последнем случае тензорный параметр разрушения фу можно заменить скалярным параметром ф кк. Далее рассматривается именно этот вариант модели.

Хотя в своем первоначальном варианте модель предложена для описания поведения горных пород, авторы работы [5] отмечают, что заложенные в ней идеи пригодны для описания поликристаллических металлов. В зависимости от выбранного для анализа структурного и масштабного уровня [1] в качестве микротрещин (понятия, с которого начинается построение модели Драгона-Мруза) можно рассматривать дефекты другой физической природы: скопления дислокаций, дислокационные петли, дислокационные субструктуры и так далее. Дефекты структуры (например дислокации) могут быть заданы геометрическими параметрами исходной (существующей до начала нагружения) ориентации и линейных размеров, а в конечном счете, локальными начальными деформациями элементарного объема. Обусловленные этими дефектами особенности процессов деформации, например поля упругих напряжений, могут быть найдены по результатам решения уравнений механики для тела с равномерным или же неоднородным распределением объемных начальных пластических деформаций.

2.1. Определяющие соотношения

Формулировка определяющих соотношений упругопластической модели Драгона-Мруза [5] базируется на следующих допущениях. Полные приращения деформаций элементарного объема (материальной точки сплошной среды) на каком-либо шаге нагружения складываются из упругих и пластических, обусловленных необратимыми сдвигами и микроразрушениями, деформаций:

dеij = dе У + dе Р. (2)

Упругие деформации связаны с напряжениями обобщенным законом Гука, причем тензор упругих постоян-

ных Бщ не изменяется в процессе микрорастрескивания:

йе і - БукІ ■ йа к

(3)

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

/(а„, к, в’ - 2 ■ Jг - 2 ■ р(к, в’■ (I0 - Jx’,

(4)

где J1 - акк, J2 - ~ — первый инвариант тен-

зора напряжений и второй инвариант девиатора напряжений соответственно; Д — постоянная материала, имеющая размерность давления; р(к, в’ — материальная функция, описывающая изменение поверхности текучести в процессе нагружения и так же имеющая размерность давления.

Таким образом, рассматриваемая модель принадлежит к классу моделей с функцией пластичности / - /(Jl, J1’, зависящей от первого и второго инвариантов тензора напряжений, и учитывающих влияние гидростатического напряжения на вид а-е-диаграммы («предел текучести») материала.

Для функции р(к, р’ в упрощенном варианте модели принимается выражение:

Р(к, в’ - Ро + к - в2 , (

где р0 — начальное значение функции.

Действительному состоянию материала {а, к, в} соответствуют условия равенства нулю функции текучести (условие текучести) и ее полного дифференциала:

/(ау, в, к ’ - 0,

й/--^— ■ йау ■ йк +-д/ ■ йв- 0, (6)

7 у Эк Эв

то есть предполагается, что, во-первых, точка {а, к, в} лежит на поверхности текучести и, во-вторых, истинное состояние материала является стационарной точкой пластического потенциала.

Ассоциированный закон пластического течения имеет вид

рі - 1, ,, (а г, в, к’ т\

йе у - Я ■ пу * '■ч * йакІ (7)

да кі

/ дау

— нормированный на-

V/дарч ■д//Эс

правляющий тензор, определяющий в пространстве напряжений внешнюю нормаль к поверхности текучести; h — функция упрочнения.

Предполагается, что степень поврежденности и разупрочнение материала нарастают пропорционально пластическому изменению объема:

рі

йв - йф кк - а ■йе р

(8)

где а — положительная постоянная для данного материала.

В свою очередь, упрочнение определяется работой напряжений и связано с девиатором пластических деформаций dejpl = dер - 3- 8у • dесоотношением:

йк - я у ■ йер1.

2.2. Упрочнение и накопление повреждений

Для проведения расчетов необходимо выразить приращения параметров поврежденности среды dp и упрочнения dк через приращения напряжений на очередном шаге нагружения, а также найти выражение функции упрочнения h.

Применяя соотношения (7), (8), для dp получим выражение:

йв - а ■ йер1 -8у - а ■ Н ■ пу ■ Э/(а, в, к’ ■ йак, ■ 8

н у у у дакІ к1

У .

(10)

Нетрудно установить, что если удовлетворяется условие пластичности (6), то имеет место соотношение:

Л да у)

пу -8у =■

-8у -

д/э//дарч ■¥/дарЧ

- 6- р(к, в’

V3//дарч -3//Зарч • (11)

Следовательно, приращение йв определяется формулой:

6- а ■ Н ■ р(к, в’ \ ,

йв- -(/Эа кі) • йа кі. (12)

-у//Эарч ■д//Эар

При выводе этих соотношений использовались следующие условия: функция р(к, в’ не зависит от а у , выражения для производных функции пластичности имеют вид

Э/(а.Г в, к’ - 2 ■ (Яу + 8у -р(к, в’’, (13)

у

а также тождество Яу ■ 8у = 0 .

Подобным же образом можно получить выражение для приращения параметра упрочнения йк . Из формул (7) и (9) следует, что

йк - ■ йер1 - ■ (Э//Эау )• Н -((/Эак )йа

у у у ^ д/Э//Эарч -Э//Эа

(14)

Сравнивая выражения (12) и (14), устанавливаем, что параметры упрочнения и накопления повреждений связаны формулой:

а рч '

Я а ■ Я а

йк -- —1— ■ йв .

3-а ■ р(к, в’ ' ' (15)

Из соотношений (4) и (6) следует, что в процессе пластического течения

Я у -*у - 2 ■ p(к, в’ ■ (10 - акк ’ . (16)

Следовательно, формула (15) сводится к выраже-

нию:

йк-2-

3а

йв .

(17)

(9) При действии отрицательного гидростатического на-

I0

пряжения (акк < 0) имеет место соотношение — <

т0 —

Ь -а кк

Для минимального упрочнения формулу а

(17) можно заменить выражением

2 і0

йк- йв . (18)

3а

Из выражений (18) и (5) для функции р(к, в’ можно получить

р(к, в’ - Р0 + —-в-в2 .

3 а

(19)

Теперь можно найти выражение для функции упрочнения h. Из второго уравнения (6) и формул (12) и (14) следует, что

Эа кі

■йа кі ■ •

д/

1 + Н-

д/ ((/ да у )у

дк д/э//да рч ■ / да

6- а ■ р(к, в’

рч

дв д/Э// Эа рч -Э// Эс

рч

-0.

(20)

Это условие должно выполняться как при нейтральном ^ - dа ы = 01, так и при активном | - da н > 0

“Л! 7------X------------ -ч

да Ы ) ^да к!

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

1

-1

Н д/Э//ЭаРЧ ' д/1Эа

х<

РЧ

э/

+дв- 6- а-Р(к ’ в’г. (21)

Как следует из формул (4), (5), на поверхности текучести должны выполняться равенства:

= -2 - (I,0 - JlУ = -2-(Ю - J,).

Эк Эк

Ж.=-2-у»-Jt).Мк_£> = 4-р-(-J,). (22)

Эр ^ 1 ,; Эр и V 1 и

Учитывая эти формулы, а также выражение (13), из формулы (21) получим:

1

10 _ т

*1 т1

Н д/Э//ЭаРЧ ■ Э//Эс

РЧ

х

х{4 • -Яу - 24 - а - в - р(к, в’} . (23)

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

н - ^ - 24-(І0, - ^’-а -р(к, в’ Н 4д/1ЭаРЧ 'д//Эарч

1( її - ^’

-в

(24)

Учитывая предположение о минимальном упрочнении материала, выражение (24) логично заменить следующим:

н - ^ - 24-(І10 - ^’-а -р(к, в’ Н -у/Э//Эарч 'д//Эарч

1 І 0 - ^-в 3а

(25)

2.5. Касательная линеаризация определяющих соотношений

Выразим приращения напряжений через приращения полных деформаций в материальной точке упрочняющейся среды. Используя формулы (2), (3), (7), можно записать следующую систему уравнений для деформаций:

йеу - V •йакі + Н - пу

Э/ (а га, в, к ’

да кі

йер1 - н-пу- д/(а”, в к’ - йак1■

і і Эан

-йа к

(26)

Умножая первое из уравнений (26) на тензор упругих постоянных Сурч, обратный тензору податливости материала 8укі, а второе — на свертку вида Сурч-прч, систему уравнений (26) можно преобразовать к системе вида:

СУРЧ - йе у - СУРЧ ' Букі -йа кі +

, 7. ^ - д/ (а ГЯ , в, к’

+ Н-Сурч ' пу----------------д---------------йа к

да кі

прч - Сурч - йеу - Н - прч " Сіірч - пу х

(27)

д/ (а га , в, к ’ да кі

-йа

кі ■

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

С- - dе(Р1) = С ^dе■■- С - dе(е1) =

^гурд ™с'у УРЧ У УРЧ у

= Сурд - dеу - Сурд - *^ук/ - dак1 . (28)

С учетом этой формулы система (27) может быть преобразована к виду:

тк/ +

Э/ (а га, в, к)

СУРЧ - йе і/ СУРЧ - Букі * йа к

+Н - Сурч - пу

Эа

--йа

кі ,

кі

п РЧ - Сіурч - йе іу п РЧ - Сіурч - Біукі - йа кі

РЧ ^УРЧ

РЧ

Н - прч - СУРЧ - пу

УРЧ Букі

д/ (а га , в, к ’ Эа кі

(29)

С учетом взаимообратности тензоров упругих постоянных Сурч и тензора податливости материала Буу, а

также определения тензора п у (формула (7)) уравнения (29) преобразуются к виду:

с йе - , н д/(ага, в, к’ с п

Сіурч * йеіу йа РЧ + Н - ~ - Сурч * пкі * йакі ,

-"УРЧ

РЧ

Эа

п - С■■ - йе ■■ -

пРЧ УРЧ У

1 + Н^п- С

1 -Г н Пу Сурч

Эу(а га , в, к ’

Эа РЧ

- «кі •йакі. (30)

Из второго уравнения этой системы можно найти выражение для скаляра пы -dаы. Подставляя это выражение в первое уравнение системы (30) и учитывая выражение (25) для функции упрочнения h, а также предполагая, что h ф 0, получим формулу касательной линеаризации определяющих соотношений для модели Драгона-Мру-за:

йа у СУРЧ - йе РЧ ,

где

"УРЧ

■С- - а

СУРЧ

эсТи ■Скірч) (эсд!„ ■Сути Iх

24 - (ї0 - Т1’- а - р(к, в’

( 1 Т 0

1 --в

3 а

, у С ^ з Стпкі

аа^

У

Эа

(31)

(32)

1

кі

а* -

0, / (ау, в’ < 0 или

/ (а у, в’ - 0,

Э/

Эа У <

1, / (ау, в’ > 0 или

/ (а у, в’ - 0,

- йау > 0.

2.4. Приведение напряжений на поверхность текучести

При правильном выборе шага по времени (или шага по параметру нагрузки в квазистатических задачах) численное решение упругопластической задачи на основе определяющих уравнений (31) не вызывает затруднений до тех пор, пока напряженно-деформированное состояние соответствует восходящей ветви а-е-диаг-раммы. В области неустойчивого деформирования, когда рост пластических деформаций продолжается при одновременном уменьшении интенсивности напряжений (нисходящая ветвь а-е-диаграммы), прямое использование формул (31) становится невозможным. Объясняется это тем, что после прохождения максимума на а-е-диаграмме функция упрочнения (25) становится отрицательной (в> !,°/3а ). Физически это означает, что способность материала к деформационному упрочнению исчерпана и превалируют процессы прогрес-

сирующего разрушения с образованием микро- или мезотрещин: начинается разупрочнение материала. На этом этапе необходимо изменение алгоритма расчета для продолжения решения. Реализованный в данной работе вариант продолжения решения основан на том, что при приближении к максимуму напряжений на а-е-диаграмме функция упрочнения стремится к нулю: Н ^ 0. Это значение фиксируется и на последующих шагах по времени не изменяется. Равное нулю упрочнение соответствует модели идеально упругопластического тела. Преобладающим процессом теперь становится накопление повреждений, которое описывается параметром в. Принимается допущение, что с этого момента все увеличение объемной деформации обусловлено необратимыми микроразрушениями. Величина dp на последующих шагах нагружения может быть рассчитана непосредственно по формуле (8). Для продолжения расчета теперь необходим алгоритм приведения напряжений на поверхность текучести (коррекция напряжений) на каждом шаге по времени (или по параметру нагрузки) для случая, когда поверхность текучести сжимается в процессе нагружения.

Для гладкой поверхности текучести эта задача решается просто. В пространстве напряжений задана поверхность текучести вида (4). Предположим, что точка а у в пространстве напряжений, изображающая напряженное состояние после очередного шага нагружения, оказалась вне области, ограниченной поверхностью текучести (4). Нормаль к этой поверхности, восстанов-

ленная из точки а у , определяется уравнением:

а у -а у д// ЭаУ |а

(33)

где X — параметр, задающий движение изображающей точки вдоль нормали, а напряжения а у принадлежат поверхности (4) и определяют точку пересечения нормали с этой поверхностью. Напряжения а у требуется найти из совместного решения уравнений (4), (33). С учетом ранее полученных формул эту систему можно представить в виде:

- 2 - Р(в) - [[ - (1 - X01 - а*кк ] = °

Яу - Яу + Х0/3 - 8у • акк

X 0 - X

2 -[у + Р(в’] 6- р(в’

■-X,

(34)

а кк

акк - акк '(1 -Х0’, где ву, ву — компоненты девиаторов напряжений в

состояниях а ¡у , а у соответственно. Из второй группы уравнений в системе (34) следует, что искомые формулы коррекции имеют вид:

- а - ,

1 1 * 1 -а

“акк - “акк _Р(в’----------,

3 3 а

(35)

а-

1

1 + 2а '

По физическим соображениям неизвестный корректирующий множитель а не должен обращаться в нуль. После подстановки этих выражений в первое из уравнений (34) получается кубическое относительно а уравнение вида:

а3 + 3-Р-а + 2-2- 0,

(36)

где

Р --2- Р(в’

І10 -акк -3-Р(в’

3-[р(в’]2 ту

Достаточно в этих уравнениях формально приравнять

, а 2 - 0, чтобы получить известные фор-

Р 2-70

3- J2У

мулы коррекции напряжений для идеально упругопластической среды с пределом текучести Y0 [8]. При численной реализации требуются аккуратный анализ возможных вариантов в зависимости от знаков коэффициентов и обоснованный выбор корней. Правильным выбором шага нагружения в расчетах практически всегда можно добиться того, чтобы разность напряжений а у - а у была мала. Кроме того, предполагая реальной оценку интервала изменения функции р(в) как 0.1 < <р(в) < 10 МПа (стократное упрочнение), можно получить соответствующую оценку коэффициентов Р, Q в уравнении (36): Р<0, Q < 0, Р < Q. При выполнении этих условий физически обоснованным значением множителя а в процедуре коррекции напряжений (35) будет первый корень уравнения (36):

а = а, = -2 - г - соз(ф/3) > 0,

/ 3 (37)

0 < соб(ф/3) = <2!г < 1,

где г -

-_/й=

2-р(в’-(І10 -акк -3-Р(в’’ 3-Ту

Анало-

гично можно рассмотреть варианты для других значений коэффициентов Р, Q. Для частного случая плоского напряженного состояния при коррекции необходимо учесть дополнительное условие а33 - 0 .

3. Методы расчета

3.1 Метод расчета квазистатических задач

Вариационно-разностный метод расчета квазистатических задач основан на вариационном уравнении

Лагранжа инкрементальной теории пластичности [6, 9]. Это уравнение имеет вид:

Ш (а ЕЕ + А*а j)^8Лej)^dV(п 1 -

-Ш( + ДР-)8(А«іУ^(п’ -

V

-¡¡( + АТ)8(АИ,-)• йS(п’ - 0 ,

(38)

где р, АРі и Т, АТі — соответственно объемные и поверхностные силы и их приращения на п шаге нагружения; Аиг- — приращения перемещений, а 8(Аиг- ’ — их вариации; а ? +А а у — модифицированный тензор напряжений Кирхгофа [9].

Модифицированные тензор деформаций Грина и тензор материального поворота имеют вид:

2 - А еіу - Аи - й,у +Аы] - й„■+

1-(ик -й, і У (ик -й,у ),

2 - А ю у = Аиг - d,у -Аи у - dн -

-(ик -d„■ ))к -d,у), (39)

где Аиг -d,у = э(Аиг )дХу , а Ху — лагранжевы координаты для (п + 1) шага нагружения. Тензор напряжений Кирхгофа аЕ +А а у преобразуется в тензор напряжений Эйлера-Коши аЕ + Аа у для (п + 1) шага нагружения по формулам:

Е -1 ЭХ/ ЭХ/ Е * \

ау + Аау = Г -^--^Х7"(а( + А ак!), (40)

г Э(х;, х2, х3 ’ йv(п+1’ тт

где Т ------1---2---— ------ . Новые координаты

Э( х 1, х 2, х 3’ йv(п’

X' - X і + Аиі, найденные на текущем шаге нагруже-

ния, становятся сопутствующими координатами для следующего шага.

В инкрементальной теории пластичности [9] используются линеаризованные определяющие соотношения вида (31), в соответствии с которыми тензор Сук/ зависит от упругих постоянных материала, вида функции пластичности и параметров, определяющих упрочнение и накопление повреждений.

Вариационно-разностная схема решения задачи (38) строится следующим образом [6, 10]. Расчетная область разбивается на шестигранные (в трехмерной задаче) или четырехугольные (в двумерных задачах) ячейки. Каждая такая ячейка рассматривается как составной элемент, собранный из тетраэдров или треугольников, рис. 1, 2. В пределах тетраэдров пространственные производные вектора приращений перемещений Аик аппроксимируются в соответствии с [8, 11, 12] по формулам:

¡¡Аик -йБ

АS

ш ^

АV

■-Аи:

-а, У-і],

(41)

где 1у — единичные орты системы координат ОХ 1Х2Х3; dS — векторный элемент площади поверхности АS, ограничивающей объем тетраэдра АV, и заданный внешним направлением нормали, а совокупность вели-

Ч Аsд

чин А,у- = ~А^~ определяет конечно-разностный аналог

дифференциального оператора d,у в пределах тетраэдра. В двумерной задаче вместо векторного элемента площади используется векторный элемент длины контура треугольника, а роль объема играет площадь поверхности, рис. 2. Индекс q введен здесь для нумерации вершин ячейки сетки, в которых вычисляются перемещения. Вычисление компонент АБд вектора пло-

щади осуществляется по формулам векторного анализа с соблюдением правила обхода контура против часовой стрелки, если смотреть со стороны внешней нормали. Так, для боковой грани с вершинами (4,1, 2) тетраэдра, заданного тройкой векторов а, Ь, с (рис. 3), имеет место формула:

AS = AS]

= 4, 1, 2

1 r г 1 = — a xb= — •

2 2

г г г

41 42 43

aX1 aX2 aX3

Ьхх bX2 bX3

(42)

Условие независимости вариаций 5(ДЦР ) перемещений в произвольном узле сетки приводит к системе конечно-разностных уравнений вида:

С* ;;

І8в

2 vp

1 ~ в

2 kp

•AVPn)-

- ( + APpp) A Vp(n) - ( + AT в ) ASPn) = 0 . Конечно-разностные операторы 5^lp, 5Цр

(43)

вычис-

ляются с учетом формул (39), (41), (42) и определяют компоненты тензоров деформаций и их вариаций через перемещения узлов сетки. Варьирование перемещений 5(Див ) в узлах сетки осуществляется с учетом граничных условий вида:

Au

X,-

= Au в (P), P(Xj ) є Su , nj = T(P), P(Xj) єS0 ,

(44)

ч X, .у, 2.

Отличие схемы аппроксимации пространственных производных, реализованной в данной работе, от схем, описанных в [8, 11, 12], заключается главным образом в том, что как в двумерных, так и в трехмерных задачах применялось двойное разбиение четырехугольной или шестигранной ячейки сетки треугольниками или тетраэдрами соответственно, рис. 1, 2. В качестве коэффициентов разрешающей системы уравнений использовались средние значения коэффициентов, полученные при первом и втором типах разбиения. На каждом шаге

нагружения при вычислении площадей треугольных элементов или объемов элементарных тетраэдров (рис. 3) осуществляется проверка знака, что позволяет контролировать ситуацию при начале «выворачивания» ячеек. Для решения системы вариационно-разностных уравнений (43) применялся метод исключения Гаусса.

3.2. Метод расчета динамических задач

Для решения задач о деформации в динамической постановке использовались уравнения сохранения массы, импульса и энергии в дифференциальной форме. В лагранжевых переменных Xj они имеют вид: pJ = const,

д v1 1 daE

' (45)

д t dW d t

p д X

j

yDij

P J ,

где р — плотность материала; J — якобиан преобразования в уравнении (40); t— время; V1 — компо-

Е

ненты вектора массовой скорости; а у — тензор напряжений Эйлера-Коши; Ь — компоненты вектора объемных сил; Ж — внутренняя энергия.

Тензор скоростей деформации Dij имеет вид:

2 • Dij = V1 • d,у + Vу • d, і. (46)

Система уравнений (45) решалась для начальных и граничных условий вида:

р(Р, 0) = ро (Р), Р(х, ^, z) є V,

Ж( Р,0) = 0, Р(х, ^, 2) є V,

a j (P,0) = 0, P(x, y, z) є V, ув (P, 0) = 0, P(x, y, z) єГ, уe (P, t) = Уe (P, t), P(x, y, z) є Su,

.E ,

(47)

Сг/ (P, г) ■ = Т (P, гX P(X, у 2) £ 50 .

В динамических задачах также были использованы определяющие соотношения в форме уравнений касательной линеаризации (31) для расчета напряженно-

a

Рис. 2. «Красно-черное» разбиение четырехугольного элемента сетки в двумерной задаче

Рис. 3. Схема вычисления площади граней и объема тетраэдра с вершинами (1, 2, 5, 4)

Рис. 4. Схема расположения узлов вторичной сетки (кружки) для

задания напряжений относительно узлов основной сетки (темные Рис. 5. Схема двухосного нагружения образца в двумерных и трех-

квадраты) мерных расчетах

деформированного состояния на восходящей (устойчивой) ветви ст-е-диаграммы и алгоритм приведения напряжений на поверхность текучести для этапа неустойчивого деформирования.

Численный метод решения динамических задач, который применялся для расчетов, подобен методам, описанным в работах [8, 11, 12]. Как и в квазистатической задаче, основные особенности разностной схемы связаны с применением аппроксимирующих соотношений для пространственных производных, построенных при двойственном разбиении расчетной области конечноразностной сеткой как для компонент тензора скоростей деформации, так и для градиентов напряжений. Для аппроксимации последних была использована вторичная сетка, к узлам которой отнесены напряжения, деформации, плотность, энергия и искусственная вязкость. Вторичная сетка сдвинута на половину шага по каждой координате относительно основной, в узлах которой определяются компоненты вектора скорости, рис. 4. Так как рассматривались случаи нагружения с малыми скоростями и небольшими давлениями, влияние искусственной вязкости было незначительно. Выбор шага по времени осуществлялся исходя из условий устойчивости разностной схемы [8, 11].

4. Численные результаты

В работе [5] приведены результаты расчета по формуле (31) ст-е-диаграммы при двухосном сжатии скальной породы с параметрами модели, приведенными в табл. 1. С целью тестирования разработанных алгоритмов были проведены расчеты в двумерных квазистати-ческой и динамической (плоское напряженное состояние), а также трехмерной динамической постановках подобной задачи. Схематически расчетный образец показан на рис. 5. Рассматривались пластинка единичной толщины размерами 1 X 1см и кубический образец с длиной ребра 1 см. Расчеты двумерных задач проводи-

лись с равномерным шагом сетки Дх = Ду = 0.04 см, а в динамической трехмерной постановке задачи — с шагом Дх = Ду = Дг = 0.1 см (количество ячеек сетки 10x10x10) и для проверки сходимости схемы с шагом Дх = Ду = Дг = 0.04 см (количество ячеек сетки 25х25х х25) при скорости деформаций |е 22| = |е 33| =4.0х101 с-1. Для однородных образцов в этом примере результаты расчета с различным пространственным шагом практически совпадают, но время расчета динамической задачи весьма значительно увеличивается при его дроблении. В квазистатических задачах рассматривался процесс последовательного нагружения с постоянным шагом перемещения противоположных граней |Ди21 = = |Ди31 = 2.0х10-5 см (рис. 5), что соответствует приращению деформаций |е221 = |езз | = 4.0 X 10-5 на каждом шаге. На следующих далее рисунках знаки компонент тензора напряжений заменены на противоположные, чтобы сохранить схему представления результатов, принятую в [5]. На рис. 6 приведены зависимости сжимающих напряжений от деформации е11. По оси ординат отложена осредненная по объему образца интенсивность напряжений стг- = (3/2^^у ■ sy) ■ ДV¿ /V , кото-\ к

рая в рассматриваемых двумерных задачах (плоское напряженное состояние) с точностью до знака совпадает с величиной сжимающих или растягивающих напряжений: стг- = |ст22| = |ст33|. В трехмерной задаче эти условия, вообще говоря, не выполняются вследствие

Таблица 1

Параметры модели для скальной породы

Плотность Модуль Коэффициент р0 , /,0, а ,

рх 10-3 , упругости Пуассона МПа МПа (МПа)1/2

кг/м3 Е, МПа V

2.7 1 471 0.3 0.095 79.33 51.42

8 =г -0.1 -0.05 0 0.05 0.1 0.15 Є,,

22 33 33

Рис. 6. Диаграмма нагружения при двухосном сжатии скальной породы: 1 — двумерный квазистатический расчет; 2 — двумерный динамический расчет; 3 — трехмерный динамический расчет; 4 — данные работы [5]

запаздывания разгрузки в направлении оси Х1. Кривая 1 соответствует результатам двумерного квазистатичес-кого расчета, кривая 2 — двумерного динамического, кривая 3 — трехмерного динамического расчета при малой скорости нагружения, кривая 4 воспроизведена по точкам, снятым с зависимости, приведенной в работе [5]. Хорошее соответствие этих зависимостей, полученных разными методами, говорит об удовлетворительной точности разработанных методов расчета. Отличие зависимости, полученной в трехмерном расчете, объясняется тем, что в объеме материала поддерживаются условия трехосного напряженного состояния. Поэтому на начальных этапах деформирования рост

величины а і =

которая про-

порциональна интенсивности напряжений сдвига, отстает от величин, соответствующих двумерной задаче плоского напряженного состояния, но затем, с началом разгрузки внутренних областей кубического образца, эта задержка компенсируется быстро прогрессирующим разрушением во всем объеме материала. На рис. 7 показана зависимость напряжений от объемной деформации. Нумерация кривых здесь та же, что и на рис. 6. Из этих рисунков видно, что деформации упругого сжатия после начала пластического течения постепенно сменяются необратимыми пластическими деформациями, сопровождающимися увеличением объема при продолжающемся действии сжимающих напряжений. Таким образом, модель описывает аномальное изменение объема материала при сжатии. На рис. 8 приведены зависимости осредненной по объему интенсивности напряжений от интенсивности деформаций сдвига, вы-

численной по формуле єі =

<%)•/ V ,

рисунке повторяют результаты, приведенные на рис. 6, кривая 3 соответствует результатам двумерного квази-статического, а кривая 4 — двумерного динамического расчета для случая двухосного растяжения с такой же по модулю скоростью, как и в предыдущем примере. Такая же нумерация кривых использована и на рис. 9, где приведены зависимости напряжений от объемной деформации для этих примеров расчета. На этих рисунках видно, что диаграмма нагружения такого материала при сжатии проходит выше диаграммы нагружения при

которые дают возможность сравнить реакцию этого материала на сжатие и растяжение. Кривые 1 и 2 на этом

Рис. 7. Зависимость напряжений сжатия от объемной деформации: 1 — двумерный квазистатический расчет; 2 — двумерный динамический расчет; 3 — трехмерный динамический расчет; 4 — данные работы [5]

= ст22= азз , МПа

1 1 1 1 1 * ■ ^ 1 1 1 1 1 1

г-'' ^ - .-*■—. - .а ...<•■■’ ■ /■ - . ■ /.• V - ] *•* ч - -У ■

1 Э 8° - 1 -

] _ Э 1,1,1 . I . I .

0.0 0.05 0.1 0.15 0.20 ^ ,

Рис. 8. Диаграммы нагружения при двухосном сжатии (кривые 1, 2) и растяжении (кривые 3, 4) скальной породы: 1,3 — двумерный квази-статический расчет; 2,4 — двумерный динамический расчет Ось абсцисс — интенсивность деформаций; ось ординат — интенсивность напряжений

растяжении. Объемная деформация при растяжении продолжает монотонно нарастать в процессе нагружения с постепенно увеличивающейся скоростью после появления пластических деформаций.

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

Рассмотрим модельную задачу о плоском напряженном состоянии в первых циклах знакопеременного квазистатического двухосного нагружения того же материала и с теми же граничными условиями, что и в предыдущем примере. Такого рода нагрузки можно создать, например циклическим изменением температуры плоского образца в условиях стесненной деформации. Пусть материал вначале сжимается из недеформированного состояния (первые 250 шагов с приростом перемещений К1 = |Ли 3 2.0 х 10 5 см на каждом шаге), и затем на-

грузка периодически меняет знак на противоположный через каждые последующие 250 шагов нагружения. При таких условиях в каждом последующем цикле нагружения предварительно созданные пластические деформации одного знака частично снимаются пластическими деформациями другого знака. Однако остаточные пластические деформации сдвига за один цикл 8р1(1),

8р1(2),... и т.д., см. рис. 10, постепенно накапливаются. Эти процессы приводят к упрочнению материала и росту предела пропорциональности (начального предела текучести, точки A, E, I, C, G на рис. 10) в каждом новом полуцикле. Вместе с тем, в материале накапливаются и повреждения (объемная пластическая деформация), что приводит к постепенному снижению модуля упрочнения (участки AB, EF, и, CD, GH). Если амплитуда изменения деформаций мала и интенсивность напряжений в первом полуцикле незначительно превышает предел пропорциональности материала в недеформиро-ванном состоянии, остаточные пластические деформации за цикл постепенно уменьшаются: 8р1(2) < 8р1(1),... и т.д. В этом случае материал может достаточно долго выдерживать циклические нагрузки без заметных разрушений, хотя процесс скрытого накопления повреждений и будет продолжаться или же может стабилизироваться, то есть перейти в область упругих деформаций. На рис. 11 показан пример расчета первых восьми циклов для случая, когда через каждые 125 шагов нагружения приращения перемещений |Лм2| = |Лм3| = = 2.0 х 10-5 см меняют знак, то есть амплитуда деформаций за один цикл уменьшена в два раза. Из этой зависимости следует, что материал постепенно приспосабливается к нагрузке и приближается к состоянию, в котором его можно назвать циклически стабильным [13]. В результате упрочнения материал переходит в область упругих деформаций.

Иное развитие событий будет наблюдаться в том случае, если в результате предыстории деформирования

Рис. 9. Диаграммы нагружения при двухосном сжатии (кривые 1, 2) и растяжении (кривые 3, 4) скальной породы: 1,3 — двумерный квази-статический расчет; 2,4 — двумерный динамический расчет. Ось ординат — интенсивность напряжений; ось абсцисс — объемная деформация

поврежденность материала близка к критической и/или амплитуда нагрузки достаточно велика. Увеличим амплитуду изменения деформаций в два раза по сравнению с первым примером расчета циклического деформирования. Тогда остаточной пластической деформации сдвига и объемной пластической деформации за первый цикл оказывается достаточно для того, чтобы при повторном нагружении в полуцикле сжатия за предел текучести (точка Е на рис. 12) перевести материал в состояние прогрессирующего разрушения (участок EF

= СТ22= а33 ’ МПа

- Полуциклы сжатия | А У F

- А 5РК1) О Л^2)/ // /о / / / / В

/ 0.02/ / е;

О 0.01 / 0.03

0 н в С Полуциклы растяжения

Рис. 10. Диаграмма двухосного циклического нагружения; упрочнение повышает предел пропорциональности (точки А, Е, I, С, О), накопление повреждений постепенно уменьшает модуль упрочнения (наклон участков АВ, ЕЕ, СО, ОН, IX)

диаграммы нагружения на рис. 12). Если предыстория деформирования материала, изображенная на рис. 12 участками диаграммы нагружения ОАВСОО, неизвестна, то переход от упругого этапа деформации (участок О Е) непосредственно к этапу неустойчивого деформирования с прогрессирующим разрушением (участок ЕЕ) будет совершенно «неожиданным». Теоретически в рассмотренной модели предысторию деформирования материала можно считать известной, если известны величины накопленной в предыдущих циклах плас-

= ст22= азз ’ МПа

Полуциклы сжатия

Полуциклы растяжения

Рис. 11. Стабилизация материала при циклической деформации малой амплитуды: материал приспосабливается к виду нагружения и переходит в область упругих деформаций

СТі = СТ22= С>зз. МПа

10

5

0

-5

-10

Рис. 12. Диаграмма двухосного циклического нагружения с большой амплитудой деформаций; упрочнение повышает начальный предел текучести для повторного сжатия (точка Е) по сравнению с пределом текучести при накоплении из недеформированного состояния (точка А), накопление повреждений приводит к прогрессирующему разрушению в начале второго цикла (кривая EF)

тической деформации сдвига 8° и объемной пластической деформации 8¥, схематически показанные на рис. 8, 9. Основное значение имеет величина 8¥, которую можно рассматривать как степень разрыхления структуры материала по сравнению с некоторой идеальной начальной структурой. Если эта величина и механические характеристики материала, подобные тем, что приведены в табл. 1, известны, то возможно адекватное описание деформационного поведения материала с учетом предыстории его нагружения как на этапе устойчивого деформирования, так и при глобальной потере сдвиговой устойчивости и прогрессирующем разрушении.

5. Заключение

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

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

На модельных примерах показано, что реализованная модель пластичности описывает различия в реакции

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

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

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

Работа выполнена при поддержке Российского фонда фундаментальных исследований, грант № 96-0100902.

Литература

1. Физическая мезомеханика и компьютерное конструирование материалов / Под ред. В.Е. Панина: В 2-х т. - Новосибирск: Наука, 1995. - Т.1.- 298 с., - Т.2.- 320 с.

2. Бриджмен П. Исследование больших пластических деформаций и разрушения. - М.: Иностранная литература, 1955.

3. Латынина Л.А. О возможных изменениях в режиме медленных движений перед землетрясением // Динамика земной коры. - М.: Наука, 1965.- С. 149-154.

4. Кафка В. Теория медленных упругопластических деформаций поликристаллических металлов с микронапряжениями как скрытыми переменными, описывающими состояние материала // Механика. Т. 7. Проблемы теории пластичности. - М.: Мир, 1976. -С. 123-147.

5. Драгон А., Мруз 3. Континуальная модель пластически хрупкого поведения скальных пород и бетона // Механика деформируемых твердых тел. Направления развития. - М.: Мир, 1983. - С. 163188.

6. Черепанов О.И., Смолин И.Ю., Стефанов Ю.П. Комбинированная вязко-упругопластическая модель среды для численного моделирования деформации и разрушения неоднородных материалов // Физическая мезомеханика. - 1998. - Т. 1.- № 2. - С. 59-72.

7. Линь Т.Г. Физическая теория пластичности // Механика. Т. 7. Проб-

лемы теории пластичности. - М.: Мир, 1976. - С. 7-68.

8. Уилкинс М.Л. Расчет упруго-пластических течений // Вычислительные методы в гидродинамике. - М.: Мир, 1967. - С. 212-263.

9. Васидзу К. Вариационные методы в теории упругости и пластичности. - М.: Мир, 1987. - 542 с.

10. Макаров П.В., Черепанов О.И., Демидов В.Н Математическая модель упругопластического деформирования мезообъема материала с ограниченным числом систем скольжения // Изв. вузов. Физика. - 1995. - Вып. 38. - № 11.- С. 26-57.

11. НохВ.Ф. СЭЛ— совместный эйлерово-лагранжев метод для расчета нестационарных двумерных задач // Вычислительные методы в гидродинамике. - М.: Мир, 1967. - С. 128-184.

12. Уилкинс М., Френч С., СоремМ. Конечно-разностная схема для решения задач, зависящих от трех пространственных координат и времени // Численные методы в механике жидкостей. - М.: Мир, 1975.- С. 115-119.

13. Писаренко ГС., Можаровский Н.С. Уравнения и краевые задачи теории пластичности и ползучести. - Киев: Наукова думка, 1981. -496 с.