CC BY
55
УДК 539.3
А.И Голованов, Л.У. Султанов Казанский государственный университет
РАСЧЕТ БОЛЬШИХ УПРУГОПЛАСТИЧЕСКИХ ДЕФОРМАЦИЙ
ТРЕХМЕРНЫХ ТЕЛ МКЭ
Abstract
The article is devoted to elasto-plastic analysis for finite deformation, large displacements and rotations. An incremental method in formulation of Lagrangian-Eulerian type is used. The Cauchy stress and objective Jaumann rate of Cauchy is represent the stressed state. The von Mises yield criterion and projection of stress to yield surface algorithm is applied.
Работа посвящена построению методики исследования напряженно-деформированного состояния (НДС) упругопластических тел с учетом больших перемещений, поворотов и конечных деформаций. Основная идея состоит в использовании метода пошагового нагружения в рамках комбинированного лагранжево
- эйлерового описания деформации среды. Напряженное состояние описывается тензором истинных напряжений Коши - Эйлера, а скорость его изменения -производной Яуманна. Условием пластического течения служит критерий пластичности Губера - Мизеса. Используется метод проецирования напряжений на поверхность текучести с итерационным уточнением уравновешенности текущего НДС. В идейном плане данная методика близка изложенным в работах [1, 2, 5-7, 9, 10 и др.] и является продолжением работ [3, 4].
Кинематика среды
В рамках лагранжевого описания деформирования среды введем в рассмотрение: ег- - орты глобальной декартовой системы координат, относительно
которой исследуется процесс деформации (рис. 1); £,] - материальные (лагранжевы) координаты точки тела, которые в конечно-элементной реализации являются
локальными координатами; г(е,1 )= х1 (к і К - радиус-вектор материальной частицы в недеформированном состоянии; и(^1) = у і к і X - радиус-вектор материальной
д у1
частицы в деформированном состоянии; J = ёй
д х]
- изменение объема в процессе
деформирования; V = Vі е і = уі е і = Я - вектор скорости.
Будем считать, что в текущий момент времени (на данном шаге нагружения) процесс деформирования представляет собой течение среды с упругопластическими свойствами. Это соответствует подходу Эйлера, широко применяемому в механике сплошных сред. В соответствии с этим кинематику деформирования тела будем определять с помощью следующих тензоров:
дvi
И =----- е і е і = кц еі е і - тензор пространственного градиента скоростей;
ду}
а = 1
2
1
ю = — 2
И + И т
И - Ит
+
ду] ду'
еіеі = е'еі - тензор деформаций скорости;
ду] ду1 )
е'еі =®у-е'еі - тензор скоростей поворота.
Определяющие соотношения
Для характеристики напряженного состояния используем тензор истинных напряжений Коши - Эйлера Ё = а 9 е іе і. Помимо скорости изменения напряжений Ё,
будем использовать индифферентную производную тензора напряжений по Яуманну Ё7, которая имеет вид
2 = Ё — ю • Ё + 2 • ю. (1)
Считая среду изотропной, для удобства разложим скорости деформаций и
напряжений на шаровую часть и девиатор, обозначив
= 3 ^іі, а 7 =1 аіі, гі'іі = —&уЛ о, а'іі =аіі — 5 у- а7,
где 5іі - символ Кронекера.
Предполагается использования аддитивного разложения полной деформации на упругую и пластическую составляющие. Физические соотношения упругого деформирования запишем как зависимость между производной Яуманна тензора напряжения и тензора деформаций скорости в виде линейного закона Гука,
Ё'7 = 20А\ Ё7 = ЪКА0, (2)
Е Е
где G = -2(1-), К = 3(1—2“), Е - модуль Юнга, д - коэффициент Пуассона.
В этом случае соотношения упругости будут полностью удовлетворять принципу
индифферентности.
Алгоритм расчета
Процесс деформирования представляется в виде последовательности равновесных состояний. Переход от предыдущего к последующему состоянию происходит путем приращения нагрузки, на каждом шаге нагружения пересчитываются геометрия и напряжения. Методика расчета состоит в разработке алгоритма
вычисления (/+1)-го состояния при известных параметрах процесса (/)-го состояния. Исходным является уравнение виртуальных мощностей в актуальной конфигурации
|[[Е : ё^О = |[|О • ёу^О + Л Р • ёу^5, (3)
О О ^п
где О - текущий объем, 5а - часть его поверхности, на которой заданы усилия, О , Р - векторы массовых и поверхностных сил. Разрешающее уравнение получим из этого уравнения мощностей путем его линеаризации, используя следующие преобразования:
іШ =:■*« = Ж
Е : дй + Е : дй + -- Е : дй 3
іО,
і ЛЛЛ 0 • д™«=ЛЛЛ|е' *.+о 3
3
• 5у іО,
і
і
ЛЛ р • шї=лл
р + р
-
3
•дуіХ,
где дй = -
•5и],
3
= 3 і0 = ё1уу .
□
о
□
□
х
В рамках классической теории течения пластическое деформирование моделируется на основе метода проецирования напряжений на поверхность текучести, которая в [8] называется процедурой «сноса напряжений на поверхность текучести» с использованием критерия Губера-Мизеса. Ввиду того, что в результате использования этого метода полученное НДС не удовлетворяет разрешающей системе уравнений, проводится итерационное уточнения НДС текущего состояния, основанное на введение в разрешающее уравнение работы «дополнительных напряжений» Е д на возможных
скоростях деформаций. Итоговое уравнение для т-й итерации на /-м шаге нагружения имеет вид
ЛЛЛ{2£г ат :д1 й' + 9К1 ат :д1 й0 +/ Е:
д1 а1 «т - «т •д1 а +д/ ат
+
□/
/ & т ] / & т / & т
^—/Е:д1 йІіО- ЛЛЛ-Г-10•§/^- ЛЛ1Р-Т~-§/^ =
I л.и /3т .и /3т
/3т
ЛЛЛ/ о 5/ ^+ЛЛ/ р 5
(4)
ЛЛЛ[/Е :5/ а-/ о •5/ уДО + ЛЛ/ Р 5 уії
ЛЛЛ/ Е
/ т /
■5 йіО.
□
□
□
х
Фигурирующие в этих уравнениях индексы текущего шага / означают, что соответствующие величины вычисляются для О/ конфигурации, т.е. уг- =гуг-, тензор
напряжений примет значение Е=/Е и аналогично для нагрузок и их скоростей.
Особо отметим слагаемое в фигурных скобках в правой части уравнения (4). По физическому смыслу это есть уравнение виртуальных мощностей (3) в /-м состоянии,
которое должно удовлетворяться для точного решения. Поскольку это уравнение на шаге нагружения в настоящей методике точно не выполняется, то его вводим в линеаризованное уравнение в виде невязки. Опыт решения нелинейных задач в шаговой постановке (в том числе и представленной в настоящей работе) свидетельствует, что наличие такого рода слагаемых в правой части соответствующего линеаризованного уравнения препятствует накоплению ошибок и не позволяет решению удалиться от действительной кривой деформирования.
Полученное уравнение является линейным относительно скорости 1 Vт , поэтому после численной дискретизации может быть получена система линейных алгебраических уравнений для соответствующих узловых значений проекций
скоростей 1у™. Так как исследуемые процессы не имеют явного динамического
характера (ускорения не учитываются), то под временем можно понимать любой монотонно возрастающий параметр, определяющий изменение нагрузки. В таком аспекте вполне уместно принять производную по времени как отношение приращения соответствующих величин, получаемых при переходе с / состояния в состояние (/+1). Например,
Теперь, в силу произвольности выбора параметра t, вполне допустимо принять
Решение задачи теории пластического течения проводится по следующей схеме. Сначала в соответствии с процедурой проецирования напряжений на поверхность
при его невыполнении, т.е. при появлении пластических деформаций, считаем, что девиатор истинных напряжений определяется как проекции девиатора «пробных» на поверхность текучести, и в этом случае
М
М
А^1. В результате получаем из (4) уравнения для приращений перемещений Аг ит .
текучести проводится упругий расчет. Вычисляем текущие скорости 1 Vт , по которым определяем тензоры 1 dт, 1 ют; с помощью соотношений (1, 2) находим скорости приращений напряжений 1Ет, шаровую часть истинных напряжений следующей конфигурации:
/+1 Е 0 =/ Е 0 +1 Е 0 М
так называемый девиатор пробных напряжений
Е'т = Е'+1Е'т Аt;
проверяем условие пластичности
Е'т+1 = е'т °Т
Далее определяем дополнительные напряжения 1Е т, фигурирующие в уравнении (4), как разность истинных и «пробных» напряжений,
Ет+1 = е'т+1 - Е'т . (5)
(5)
По достижении сходимости, т.е при достаточно малых (5), переходим к следующему шагу нагружения. Определяем конфигурацию 1+1 шага 1+1 у1 = 'у1 +А 1п1 и соответствующее напряженное состояние 1+1 £ = £ т .
Конечно-элементная дискретизация
В качестве базового возьмем произвольный восьмиузловой конечный элемент, изображенный на рис. 2 в декартовой и локальной (безразмерной, отнесенной к размерам элемента) системах координат.
Введем аппроксимацию геометрии и скорости следующим образом:
у (' )=£'у(1)
,=1
'V (' )=£ (])
,=1
Рис. 2. Восьмиузловой конечный элемент в декартовой (а) и локальной (б) системах координат
1у\
где N (] )=1 ( + £,,V )1 + £,,2£,2 Л + £,,3£,3) - функция формы, -1 <£\ £2, £3 < 1 8
координаты узлов (' - номер шага нагружения, 1 - номер проекции, , - номер узла в элементе), 'у! - скорость узла.
Приведенные выше функции формы обладают свойствами полноты и конформности, т.е. образуют полные системы полиномов и при стыковке двух элементов обеспечивают непрерывность функций 'у , 'я..
тт д' у д' у Я б
Для перехода от —- к строим матрицу Якоби
дГ ду
Вычисляем матрицу, обратную \3 ], с помощью которой определяются искомые производные
с
1А~1 ]
д =с д ^ Сі] д? ]
д у
1л I
определяющие компоненты тензоров п , а , ю .
В результате обычных для метода конечных элементов процедур получим на т-й итерации 1-го шага нагружения систему линейных алгебраических уравнений
2
е
2
1
[1К]{а1дт }= {і Р}- }-{^т }.
Ее решение дает значения приращений перемещений А1 uт, по которым вычисляются параметры НДС по вышеописанной схеме.
Числовые примеры.
1) Рассмотрим задачу об изгибе упругой балки прямоугольного поперечного сечения в кольцо. Балка жестко защемлена с одной стороны, к другой приложен изгибающий момент. Длина балки 1=50 см, высота И=1 см, ширина 6=0,125 см, модуль упругости £=1000 кг/см , коэффициент Пуассона ц=0. Значение момента, при котором балка изгибается в кольцо, было взято из аналитического решения М=2п кг см.
Используется сетка вышеописанных конечных элементов размером 200*1x1. Нагрузка была разбита на 1000 шагов. Параметрическое исследование по изменению сетки конечных элементов показало, что увеличение числа элементов по высоте мало сказывается на точности решения, наибольшее влияние оказывает число элементов по длине. Величина шага нагружения существенно влияет на точность, что вполне объяснимо. В качестве примера на рис. 3 изображены деформированное состояние балки и несколько промежуточных этапов нагружения.
Рис. 3. Процесс деформирования балки
2) Следующей была задача об упругопластическом деформировании жестко защемленной с обоих концов балки под действием распределенной нагрузки q=21 кг/см2. Длина балки /=25 см, высота И=1 см, ширина 6=0,125 см, модуль упругости £=20000 кг/см, коэффициент Пуассона М=0, предел текучести от=150 кг/см2.
Поскольку задача является симметричной, то достаточно рассмотреть половину балки, введя дополнительные условия. При решении используем сетку конечных элементов размером 100x10x1. Нагрузку разбили на 100 шагов. На рис. 4 показано распределение интенсивности напряжений на промежуточных и конечном этапах нагружения.
Интересно отметить в этой задаче наличие впадины в окрестности защемления, а именно в этой зоне возникают конечные пластические деформации. Возможно, этот эффект имеет общую природу с появлением шейки в растянутом цилиндрическом образце из пластичной стали.
Для контроля задача была решена на известном программном комплексе АКБУБ. Результаты этих расчетов качественно и количественно совпали с полученными в настоящей работе.
Рис. 4. Распределение интенсивности напряжений в балке на промежуточных и конечном этапах нагружения
Библиографический список
1. Голованов А.И., Бережной Д.В. Метод конечных элементов в механике деформируемых твердых тел. - Казань: ДАС, 2001. - 301 с.
2. Голованов А.И., Кузнецов С.А., Расчет больших деформаций неупругих тел в комбинированной лагранжево-эйлеровой постановке // Модели механики сплошной среды: Обзорные доклады и лекции XVI сессии Международной школы. - Казань: Изд-во Казанского математического общества, 2002. - С. 5-21.
3. Голованов А.И., Сенникова А.С., Султанов Л.У. Постановка задачи и разработка алгоритма исследования больших деформаций упругих тел // Зимняя школа по механике сплошных сред (тринадцатая): Тезисы докладов. - Пермь, 2003. - С 109.
4. Голованов А. И, Султанов Л. У. Пошаговое исследование больших деформаций упругих тел МКЭ. Тезисы докладов Двенадцатой международной конференции по вычислительной механике и современным прикладным программным системам. -М.: Изд-во МАИ, 2003. - Т. 1. - С. 191-198.
5. Кукуджанов В.Н., Кондауров В.И. Численное решение неодномерных задач динамики твердого деформируемого тела // Проблемы динамики упругопластических сред. - М.: Мир, 1915. - С. 39-84.
6. Коробейников С.Н. Нелинейное деформирование твердых тел. - Новосибирск, 2000.
- 262 с.
7. Поздеев А.А., Трусов П.В., Няшин Ю.И. Большие упругопластические деформации: теория, алгоритм, приложения. - М.: Наука, 1986. - 232 с.
8. Уилкинс М.Л. Расчет упругопластических течений // Вычислительные методы в гидродинамике. - М.: Мир, 1967. - С. 212-263.
9. Taylor L.M., Becher E.B. Some computational asrect of large deformation, rate-dependent plasticity problems // Comput. Meth. Appl. Mech. Eng. - 1983. - Vol. 41. - № 3. - P. 251-277.
10. Yamada Y. Nonlinear matrices, their implication and applications in inelastic large deformation analysis // Comput. Meth. Appl. Mech. Eng. - 1982. - Vol. 33. - № 1-3. - P. 417-437.
Получено 10.10.2003
