CC BY
70
УДК 539.3
И.Ю. Смолин1, 2, Р.А. Бакеев1, П.В. Макаров1, 2
1 Институт физики прочности и материаловедения СО РАН (г. Томск) 2Томский государственный университет
ЧИСЛЕННОЕ РЕШЕНИЕ НЕКОТОРЫХ ДВУМЕРНЫХ ЗАДАЧ ДЛЯ УПРУГОПЛАСТИЧЕСКОЙ МИКРОПОЛЯРНОЙ СРЕДЫ
Abstract
Numerical solution of the plane strain problems for elastic plastic Cosserat medium is considered. Some features of the finite difference method application for the equation set governing elastic plastic deformation in the media with independent rotations are discussed. To define the time step the speeds of elastic waves connected with internal rotations are taken into account. A test calculation of plastic strain localization in homogeneous media with a stress concentrator is solved and discussed.
Как известно из экспериментов, в поликристаллах наряду с трансляционными модами деформации развиваются и поворотные моды. Вращению подвержены как отдельные зерна, так и их фрагменты, а также конгломераты зерен. К основным причинам возникновения внутренних поворотов можно отнести недостаток активных систем скольжения, взаимодействие потоков дефектов на границах фрагментов мезоструктуры. Недостаточная аккомодация сдвига в сопряженной системе скольжения в условиях стесненной деформации приводит к возникновению крутящего момента. При этом действие внутренних моментов, оказывающих существенное влияние на поведение фрагментов структуры на мезоуровне, на макроуровне усредняется и становится равным нулю. Поэтому вводить в рассмотрение моменты имеет смысл только на определенном масштабе рассмотрения, связанном с характерными размерами внутренней структуры материала. Модель микрополярной среды рассматривается именно как модель мезоуровня. Имеющаяся в такой среде дополнительная степень свободы поворотного типа является откликом неоднородного развития деформации на микроуровне. При осредне-
нии всех действующих в представительном мезообъеме силовых и моментных напряжений должны получиться напряжения, характерные для макрочастицы. Это условие определяет связь мезоуровня с макроуровнем.
В работе предложен способ численной реализации уравнений, описывающих движение среды Коссера. Рассмотрены особенности численной реализации модели: центрирование разностной схемы, выбор расчетного шага по времени, ограничение моментных напряжений. Приведены примеры расчетов упругопластического поведения микрополярной среды, полученные результаты сравниваются с решениями для классической модели.
Математическая постановка
Запишем полную систему уравнений моментной модели среды Коссера [1-3] для случая плоского деформированного состояния.
Закон сохранения количества движения для случая отсутствия массовых сил имеет обычный вид
Закон сохранения моментов количества движения при отсутствии массовых распределенных моментов записывается следующим образом:
где Оу - компоненты тензора силовых напряжений, щ - компоненты тензора моментных напряжений, иі - компоненты вектора перемещения, шг- - компоненты вектора независимых поворотов, р - плотность материала, у - плотность момента инерции частицы при вращении.
В среде с поворотами тензор полной деформации несимметричен у у Ф у и выражается через симметричный тензор деформации
еіу и повороты ю3 и й3 = - й12 = Й21 = 1/2 (й2 1 - и12) следующим образом [1]:
Тензор изгибов-кручений определяется через вектор поворотов,
°11,1 + °21,2 = P U1-> °12,1 + °22,2 = P U2 .
(1)
°12 °21 + М-13,1 + M-23,2 = j Ю3 .
(2)
Y11 = U1,1 , 122 = U2,2 , У12 = U2,1 ®3 , 121 = U1,2 + ®3
(3)
(4)
Законы сохранения массы и баланса энергии имеют вид
р0Г0 = pV, (5)
р Е = Опуп + 012У12 + °21У21 + О22у22 + М13к13 + М23к23 , (6)
где р0, р, У0, V - начальное и текущее значения плотности материала и объёма некоторой малой области материала соответственно, Е -удельная внутренняя энергия единицы массы.
Чтобы описать процесс упругопластической деформации, необходимо определяющие соотношения записать в дифференциальной форме. Это позволяет распространить подход на случай конечных деформаций и воспользоваться линейными определяющими уравнениями, которые необходимо записать в релаксационной форме. Пусть полные тензоры скоростей деформации уТ и изгибов-кручений кТ
есть суммы их упругих и пластических составляющих,
= У у = У^ + Ур кТ = к и = к е + к р (7)
Теперь полные у Т и к Т вычисляются по формулам (3) - (4), а их
пластические составляющие следует задать любым удобным способом, вытекающим из рассмотрения процесса пластической деформации.
С учетом (7) определяющие уравнения могут быть записаны для изотропного центрально симметричного тела в виде
О11 = ^ ( у кк- укк)+2 м (Уп - Уп),
О22 = ^ (Укк - Укк ) + 2 М ( У22 - У22 ),
О12 = (М + а) (У12 - Ур2 ) + (М- а) (У21 - У21 ) ,
О21 =( М+ а)(У 21 - У 21 ) + ( М- а)( У12 - Ур2 ) , (8)
М13 =( V + £ )( к1Т3- к1Р3) ,
М23 =( V + £ )( к Т3 - к 23 ) ,
М31 =( V- £)(кТ3- кр3) ,
М32 = ( V- £ ) ( к 23 - к 23 ) .
Здесь X и ц - коэффициенты Ламе; V, в - новые упругие модули; - означает коротационную производную Яумана, +
+sikО.. Компоненты тензора моментных напряжений ц31 и ц32 не
равны нулю, но не входят ни в уравнения движения, ни в уравнения баланса энергии, поскольку так же, как и компонента тензора силовых напряжений о33 , обеспечивают условия плоского деформирования.
Следовательно, в данной постановке присутствуют лишь два новых физических параметра среды а и V + в.
Пластические составляющие деформаций и изгибов-кручений для двухмерного случая могут быть представлены в виде [3]
где функции пластичности определены следующим образом:
Здесь параметр I - характерный размер микрополярной среды. С физическими параметрами он связан следующим соотношением:
V + 8 = 2 .
Выписанную систему уравнений следует дополнить соответствующими начальными и граничными условиями.
Особенности численной реализации
Задача о моделировании упругопластического деформирования среды, описываемая записанной выше системой гиперболических уравнений, в общем случае не имеет аналитического решения. Для ее численного решения в работе используется способ построения разностной схемы, основанный на представлении частных производных по пространственным переменным контурными интегралами, применении теоремы о среднем интегрального исчисления и аппроксимации производных по времени обычным конечно-разностным методом. Этот метод известен как метод Уилкинса [4]. Получаемая разностная схема является схемой второго порядка точности.
13 М 0Ц ’ 23 М Ям ’
13 23
/ (му-) = >/1,5(М123 +м23) - ¥ 1
При решении задачи для классической симметричной среды компоненты координат и векторов скоростей определяются в узлах расчетной сетки, а плотность, внутренняя энергия, компоненты тензоров напряжений и скоростей деформации относятся к центрам ячеек. В схеме Уилкинса аппроксимация производных по пространственным переменным для функций, определенных в узлах и ячейках расчетной сетки, строится по-разному. Если функция определена в центрах ячеек, то ее частная производная по пространственным координатам будет определяться в узлах расчетной сетки. И, наоборот, для функций, определенных в узлах сетки, частная производная по пространственным координатам будет определяться в ячейках.
Для построения конечно-разностных уравнений модели среды Коссера за основу была взята уже готовая схема для классической среды. Ее необходимо было дополнить конечно-разностными аналогами геометрических соотношений для компонент тензора изгибов-кручений, определяющих соотношений для моментных напряжений и уравнения закона сохранения момента количества движения. Кроме того, были модернизированы соотношения для недиагональных компонент тензоров силовых напряжений и деформаций, а также закона сохранения количества движения с учетом несимметричности тензоров деформаций и силовых напряжений.
В зависимости от того, к узлам сетки или к центру ячеек будут отнесены новые переменные (компоненты поворотов и моментные напряжения), существуют два способа построения конечно-разностных уравнений. Это обусловлено видом закона сохранения момента количества движения, в левую часть которого входят силовые напряжения, определенные в ячейках, и частные производные по пространственным координатам от моментных напряжений. Необходимо отцентрировать это уравнение, т. е. привести его либо к узлам сетки, либо к центрам ячеек. Формально необходимо решить вопрос: в узлах или в ячейках определяются независимые повороты.
Первый способ: повороты определены в ячейках, моментные напряжении - в узлах. В этом случае уравнение центрируется в ячейках и исчезает необходимость дополнительных осреднений. Второй способ: повороты определены в узлах, моментные напряжения - в ячейках. Центрирование происходит в узлах сетки, значение силовых напряжений в узлах определяется как среднее по окружающим ячейкам. Авторами были реализованы и проанализированы оба способа аппроксимации и согласно физическим соображениям вы-
бран второй. Это обусловлено тем, что в случае наличия классической среды уравнение движения центрируется в узлах, т. е. сплошная среда как бы заменяется набором дискретных масс, сосредоточенных в узлах. Векторы скорости приписаны к узлам. Логично предположить, что эти массы могут и поворачиваться, тогда вектор поворота тоже следует приписать к узлам расчетной сетки.
Отдельного рассмотрения заслуживает условие устойчивости полученной разностной схемы. Классическое условие устойчивости, предложенное фон Нейманом и Рихтмайером, имеет вид
где с - скорость звука, а Аt и Ах - шаги по времени и пространству соответственно. Предложено несколько обобщений этого условия для расчета, например, ударных волн или многомерных задач, но они касались только классической среды. Наибольшую скорость в такой среде имеет продольная упругая волна, и именно эта скорость используется при условии устойчивости. Однако в среде Кос-сера возможно существование новых типов упругих волн, которые отсутствуют в классическом случае. Скорости этих волн могут превышать скорость продольной волны, поэтому их необходимо учитывать для получения устойчивого решения.
Как показано Новацким [1], в упругой изотропной среде Коссе-ра может существовать 4 типа упругих волн. В зависимости от условий задачи необязательно все они будут реализовываться одновременно. Только одна из этих волн не является диспергирующей. Это продольная волна, которая существует и в классическом случае. Ее фазовая скорость, как и в симметричной теории упругости, задается выражением
с-=^.
Фазовая скорость волны кручения (или волны микровращений) задается выражением
с— < 1, Дх
с
где ш
плотность момента инерции
элементарного объема среды, а г - радиус инерции.
Этот тип волн не существует в классической среде. Ее скорость зависит от частоты ю, т.е. волна этого типа обладает дисперсией. Кроме того, волна микровращений может существовать только при достаточно больших частотах ш > ш *, и в пределе очень больших частот (при Ш —— да ) ее скорость стремится к С3 .
Для определения скоростей сдвиговых волн нужно проанализировать соответствующее дисперсионное уравнение [1]
Это дисперсионное уравнение является биквадратным уравнением относительно ю (или к). Аналитическое решение этого уравнения удается записать, но получить выражения для фазовых скоро-
сложно. Однако можно построить численное решение этого уравнения и проанализировать его.
Дисперсионное уравнение (9) определяет две дисперсионные ветви в плоскости ю - к, которые характеризуют две сдвиговые волны. Обе эти волны являются диспергирующими. Характерной особенностью одной из кривых является то, что она не выходит из начала координат и соответствующий ей тип сдвиговых волн не может распространяться с частотой ниже ш *, как и волна микровращений.
Расположение дисперсионных кривых определяется соотношением параметров модели и для двух характерных наборов этих констант представлено на рис. 1. Там же показаны асимптоты, определяемые тремя характерными фазовыми скоростями с2, с4, с5. Отличием этих наборов физических параметров является то, что характерная фазовая скорость с4 оказывается либо больше, либо меньше скоростей с2 и с5, а дисперсионные кривые, соответственно, имеют разные асимптоты при к — да. Поскольку для скорости с4 можно записать
ш4 -ш2[к2(с22 + с42)+ ш2]+ш*2с52к2 + с^с^к4 = 0,
(9)
йш
йк
довольно
с Л =.
у.2
рг
2
1
= с*-.
где 1 - параметр в определяющих соотношениях, характеризующий «моментную» жесткость среды, а г - радиус инерции, характеризующий её инерционные свойства, то рассмотренные случаи характеризуются тем, какая из двух характерных длин микрополярной среды 1 или г является больше.
600
500
400
300
200
100
СО
N •’З-\
' сь
200
400
600
800
Ш
400
350
300
250
200
150
100
50
ю Ш \° \
' ф /' __ / '
. X/ . ^
к \ ^ С4
у'
200
400
600
800
Рис. 1. Дисперсионные кривые для сдвиговых волн в среде Коссера и асимптоты характерных фазовых скоростей для разных соотношений характеристик среды: I > г (а), I < г (б)
Графики зависимости фазовых скоростей всех упругих волн от частоты представлены на рис. 2 для тех же двух наборов упругих констант, что и на рис. 1. Только одна скорость продольной волны
г
не зависит от частоты. Ее график - прямая линия, параллельная оси абсцисс. Остальные скорости зависят от частоты. Обращает на себя внимание особенность двух волн, скорости которых стремятся к бесконечности при приближении частоты к предельной величине ш *.
Рис. 2. Зависимости скоростей упругих волн в среде Коссера для разных соотношений характеристик среды: I > г (а), I < г (б);
1 - продольная волна, 2, 4 - сдвиговые волны, 3 - волна кручения
Анализируя расположение дисперсионных кривых для разных наборов свойств материала, можно сделать два важных вывода, которые надо учитывать при численном решении динамических задач:
1) При определенном сочетании значений упругих констант и плотности момента инерции скорость волны кручений и одной из сдвиговых волн может превысить скорость продольной волны даже при очень больших частотах ш ^ да .
2) Обе эти волны имеют неограниченно большие скорости распространения при частотах, стремящихся справа к ш *.
Поскольку при распространении произвольного импульса присутствуют гармоники разных частот, то при численном решении динамических задач в среде Коссера возможен неустойчивый счет при сколь угодно малых конечных шагах по времени (когда скорость волн стремится к да). Для соблюдения условия устойчивости следует использовать максимальную скорость из набора {с1, с3, с4}, причем для последних двух использовать еще и дополнительный множитель = 1,5, чтобы сузить область частот, в которых может наблюдаться неустойчивость.
Если две величины 1 и г, имеющие размерность длины, близки, то скорости двух характерных сдвиговых волн с4 и с5 также примерно равны и оказываются меньше классической сдвиговой волны с2. В данной статье мы будем рассматривать именно такой случай. Если при этом Сз < С1 , то условие устойчивости можно оставить прежним, как в классическом случае.
Результаты расчетов и их обсуждение
Выполнено численное моделирование одноосного растяжения однородного образца из алюминиевого сплава. Целесообразность использования однородного материала оправдана тем, что расчет носит тестовый характер и направлен на выяснение особенностей модели. В силу симметрии задачи расчет проводился для одной четвертой части образца. На верхней грани прямоугольной расчетной области задавались скорости, соответствующие растяжению, на правой грани - условия отсутствия напряжений (свободная поверхность), а на оставшихся двух гранях - условия симметрии.
Как с практической, так и с теоретической точки зрения интересным явлением при пластическом деформировании является образование полос локализованной деформации. Их появление в расчетах обусловлено свойствами определяющих соотношений и провоцируется существованием в среде концентраторов напряжений. Простейшей моделью, допускающей локализацию деформации и образование полос сдвига, является идеальная пластичность. В качестве концентратора напряжений можно задать особый вид граничных условий на верхней грани, соответствующий растяжению вдоль верти-
кальной оси и отсутствию перемещений вдоль горизонтальной оси («жесткий захват»).
На рис. 3 приведены картины распределения интенсивности пластических деформаций, полученные в расчетах для классической среды. В первом случае расчет проводился с использованием модели идеальной пластичности. Образовалась шейка с явно выраженными полосами сдвига, направленными под углом 45° к оси растяжения. Во втором случае использовалась модель с упрочнением. Применялся линейный закон упрочнения вида
Г = Г0(1 + аер1),
где ер - интенсивность пластических деформаций, а - некоторый коэффициент.
Как видно из рисунка, полосы зарождаются на жестком захвате. Затем в случае идеальной пластичности мощная полоса локализованной деформации образуется в области шейки, а в случае среды с упрочнением более деформированная область расположена в углу около «захвата».
а б в
Рис. 3. Распределение интенсивности пластических деформаций для разных моделей среды: а) классическая модель с идеальной пластичностью; б) классическая модель с линейным упрочнением; в) моментная модель с идеальной пластичностью
Результат, приведенный на рис. 3, в, получен в расчете с использованием модели среды Коссера, но также с условием идеальной пластичности. Видно, что качественно картина распределения пластической деформации схожа с расчетом по классической модели с упрочнением. Если рассмотреть количественные характеристики, то для среды Коссера характерно снижение значений силовых напряжений и накопленных сдвигов. Действие концентраторов напряжений как бы «смягчается». Таким образом, развитие независимых поворотов приводит к такому же эффекту, как и деформационное упрочнение в классическом случае. Аналогичный вывод был сделан нами ранее на основе анализа осредненных а-е диаграмм для среды Коссера [5].
На рис. 4 приведены распределения интенсивностей пластических деформаций вблизи мощного концентратора напряжений -«жесткого захвата» для образца других размеров. На рис. 4, а показано поведение классической идеально пластичной среды. Наблюдается сильная локализация сдвиговой деформации в полосе. Искажение сетки происходит вблизи концентратора и в полосе. Вместе с тем в образце можно выделить малодеформированные области, которые находятся еще на упругой стадии деформационного процесса, т.е. деформация объема носит ярковыраженный неоднородный характер. На рис. 4, б представлен расчет, выполненный для микрополярной модели. Полоса локализации в этом случае не формируется, накопление пластической деформации происходит в малой области вблизи концентратора. В остальной части образца пластические деформации распределены почти однородно. Приведенный расчет демонстрирует важное свойство модели среды Коссера. Его можно
а б
Рис. 4. Формирование полосы локализации пластической деформации вблизи концентратора напряжений в случае идеальной пластичности: а) классическая модель; б) микрополярная модель
трактовать так, что наличие дополнительных степеней свободы позволяет перераспределять энергию внутри образца. Диссипация энергии при пластическом деформировании возможна не только за счет пластических сдвигов, но и за счет независимых поворотов. Выбором значений параметров модели мы можем регулировать способность материала к внутреннему вращению. Другими словами, можно задать такие физико-механические параметры материала и условия нагружения, при которых элементарным частицам среды будет энергетически выгоднее вращаться.
Заключение
В последнее время наблюдается рост интереса к применению моделей обобщенных сред, а также к исследованиям процессов деформации и разрушения материалов на мезоуровне. Это обусловлено двумя факторами. Во-первых, бурный рост и доступность вычислительной техники стимулировали возможности численного решения задач механики, в том числе и с применением неклассических моделей сред. Во-вторых, пришло осознание того, что без учета неоднородностей внутреннего строения материалов невозможны объяснение особенностей их механического поведения и разработка материалов нового поколения. Поскольку применение микрополярных моделей сред позволяет неявно учесть микроструктуру материала, то, по сути, эти модели являются моделями промежуточного уровня между макроскопическим осредненным подходом и микроскопическим уровнем осредненного рассмотрения.
Метод численного решения, выбранный в данной работе, хорошо зарекомендовал себя как для решения динамических задач распространения упругих и ударных упруго-пластических волн, так и задач квази-статического нагружения. Применение его к модели среды с дополнительной степенью свободы поворотного типа и моментными напряжениями потребовало решения ряда задач, связанных с центрированием разностной схемы и выполнением условия устойчивости. Эти проблемы обсуждены в работе и даны определенные рекомендации.
Получены решения для однородной микрополярной среды и проанализированы отличия от соответствующего решения для классической симметричной среды. На их основе можно сделать вывод о том, что наличие в среде дополнительной степени свободы
оказывает значительное влияние на картину распределения пластической деформации.
Работа выполнена при частичном финансировании по Президентскому гранту поддержки ведущих научных школ России НШ-394.2006.1 (школа академика В.Е. Панина), а также по гранту РФИИ № 07-05-00274-а.
Библиографический список
1. Новацкий В. Теория упругости / В. Новацкий. - М.: Мир, 1975. - 872 с.
2. Эринген А.К. Теория микрополярной упругости / А.К. Эрин-ген // Разрушение. - М.: Мир, 1975. - Т. 2. - С. 646-751.
3. Forest S. Elastoviscoplastic constitutive frameworks for generalized continua / S. Forest, R. Sievert // Acta Mech. - 2003. - Vol. 160. -P. 71-111.
4. Уилкинс М.Л. Расчет упруго-пластических течений / М. Л. Уилкинс // Вычислительные методы в гидродинамике / под ред. Б. Олдера, С. Фернбаха, М. Ротенберга. - М.: Мир, 1967. - С. 212-263.
5. Смолин И.Ю. Обобщенная модель упруго-пластической среды с независимыми пластическими поворотами / И.Ю. Смолин, П.В. Макаров, Р.А. Бакеев // Физическая мезомеханика. - 2004. -Т. 7. - Спец. выпуск. Ч. 1. - С. 89-92.
Получено 29.06.2007
