Конечно-элементный решатель упругопластических задач Текст научной статьи по специальности «Строительство и архитектура»

МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ И КОМПЬЮТЕРНЫЕ ТЕХНОЛОГИИ

УДК 624.44:539.376

конечно-элементный решатель упругопластических задач

© 2007 г. П.П. Гайджуров

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

- адаптивный спуск (автоматический выбор величины шага нагружения - метод корректирующих дуг);

- процедура с обновлением матрицы жесткости (МЖ) ансамбля конечных элементов (КЭ) на каждой итерации;

- процедура обновления МЖ ансамбля КЭ на каждом шаге нагружения (внутри итерационного цикла МЖ постоянна);

- процедура с использованием начальной МЖ, не изменяющейся на протяжении всего расчета.

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

В статье [1] приведен стартовый итерационный алгоритм, построенный на основе МКЭ в перемещениях в сочетании с методом переменных параметров упругости (МППУ). Выбор МППУ главным образом объясняется принятым способом решения системы линейных алгебраических уравнений (СЛАУ) методом сопряженных градиентов с предварительным нормированием МЖ ансамбля КЭ. В публикации [1] предлагалось судить о переходе материала в пластическое состояние по достигнутому уровню интенсивности деформаций в центрах КЭ. Такое допущение объяснялось малыми размерами элементов и экономией вычислительных ресурсов.

В настоящей работе представлен вариант нелинейного решателя, базирующийся на более точном моделировании упругопластических свойств материала путем анализа деформаций в узлах интегрирования (точках Гаусса) КЭ. Напомним, что для «линейки» объемных изопараметрических КЭ с полилинейной аппроксимацией перемещений таких точек восемь. Кроме этого в предлагаемой интерпретации МППУ делается упор не на автоматический [1], а на авторизованный подход к анализу сценария нагружения и

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

Практическая реализация МППУ предполагает задание диаграммы деформирования материала ст I ~ е ^ (интенсивность напряжений ~ интенсивность

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

- И. А. Биргера [2]:

е ■

Ec =-

1 _ 1 - 2 v

2 3 E

а

е ■

1 +

1 - 2 v

3E

v c =-

а

е ■

G c =-

1+

1 - 2 v

3E

а

е

3 е

(1)

- П. А. Лукаша [3]: а ■

Ec =-

е

v c = 2

Er

1--L (1 _ 2v )

E

Gc =-

E

_e E"

(2)

-1 + 2v

с, О с, V с - переменные модуль упруго-

Здесь Е

сти, коэффициент Пуассона, модуль сдвига; Е и V -модуль упругости и коэффициент Пуассона, соответствующие отрезку диаграммы деформирования е I < е I (е ^ - величина интенсивности деформаций равная условному пределу пропорциональности). Если пользоваться схематизированной (кусочно-линейной двухзвенной) диаграммой ст ^ ~ е ^ с линейным упрочнением, то выражение для «секущего» модуля преобразуется к виду

а

Ec = E[1 -ю(е г )]е г, где функция пластичности А. А. Ильюшина [4]

О, при е i < е 1t ,

Е

ю(е i) =

f

1—

(1 -п), при еi >е, ;

Ek E П =- - параметр упрочнения; E k - «касатель-

E

ИЬШ» модуль упругости.

Соотношения (1) выведены исходя из гипотезы пропорциональности девиаторов напряжений и деформаций, а также с учетом того, что при одноосном напряженном состоянии (а { = а) интенсивность

деформаций определяется по формуле

1 - 2v Е : = Е--а.

г 3E

Выражения (2) получены в предположении, что производная удельной энергии деформации по третьему инварианту тензора деформаций равна нулю.

Общим предположением при выводе зависимостей (1) и (2) является гипотеза о выполнении объемного закона Гука.

Ес • 10-5, МПа 2,0 ч

0 0,002

0,006

0,010

0,014 е,

а)

Vc

0,46.

0,42.

0,38-

0,34-|

0,30

002 '0,0060,010 0,014 0,018 е,

0 0,002 ' 0,006 б)

"0,01 0,014 Е,

Рис. 1. Графики Ес ~ е t (а) и v c ~ е t (б): 1 - модель И.А. Биргера; 2 - модель П.А. Лукаша

Сравнение значений Ec и vc, полученных по формулам (1) и (2) для обычной стали (E = 2-105 МПа, v = 0,29, п = 0,03425, е i = 0,002), показано на рис. 1.

Как и следовало ожидать, с ростом интенсивности деформаций е i величина E c уменьшается, а коэффициент Пуассона v с увеличивается, стремясь к

значению 0,5, т. е. переходя в пластичное состояние материал, становится практически несжимаемым. Заметную разницу в результатах формулы (1) и (2) имеют на участке 0 <ei < 0,01. Это объясняется тем, что при е i = 0 переменные параметры упругости E c, Gc, vc в соотношениях (1) автоматически не переходят в параметры упругости E, G, v. Учитывая, что выражения, предложенные П.А. Лукашом, более практичны в плане построения вычислений дальнейшие рассуждения будем строить на их основе.

Основные нюансы предлагаемого алгоритма решения упругопластических задач в рамках деформационной теории отражены в блок-схеме на рис. 2. Процесс активного нагружения (nagr = 1) организован в виде последовательности шагов step = 1,2, .... При этом возможны два подхода к управлению внешними силами: аддитивное увеличение нагрузки Pk+1 = Pk + APk+1, k = 1,2,... ; дискретное приложение нагрузки «порциями» AP k , здесь индекс k соответствует номеру шага step нагружения. В частном случае AP = const для первого варианта нагружения имеем Pk = kAP .

При аддитивном нагружении результирующее поле перемещений (массив w(ng)) определяется непосредственно после успешного завершения итерационного цикла. При дискретном нагружении для вычисления результирующих узловых перемещений необходимо использовать принцип суперпозиции на серии шагов k = 1,2,.... В обоих вариантах механические характеристики материала в случае пластического деформирования на шаге нагружения k непрерывно изменяются в процессе итерационного уточнения решения и «наследуются» при переходе к k +1 шагу нагружения. Момент начального перехода материала в пластическое состояние оценивается по заданной величине е i . Для учета разгрузки следует назначить

признак задания nagr = -1.

Рассмотрим тестовые примеры, решенные с помощью исследовательского комплекса POLYGON, разработанного автором, и коммерческих комплексов ABAQUS и COSMOS.

Пример 1. Центральное сжатие куба с размерами 0,1x0,1x0,1м равномерно распределенной силой интенсивностью q (свободная осадка), см. рис. 3. Материал сталь 40Х (термообработка - отжиг) и арматурная сталь A-III, схематизированные диаграммы деформирования показаны на рис. 4.

©-Ч

Решение СЛАУ

Ввод признаков задания _step и nagr_

[H] te Г {W} = {Р}

step

Вычисление элементов массива eia (ne,8)

Нет

Ввод топологии и геометрии объекта

Ввод механических характеристик и параметров нелинейных моделей

I ~

Инициализация массивов е1а (пе,8) и е1ш (пе,8)

г

Чтение с диска массивов е1а (пе,8) и е1ш (пе,8)

Инициализация массивов w (ng) и wa (ng)

\iter = 0

^ Да ^

Корректировка массива

eia (i, j) = eiw (i, j) + nagr• eia (i, j), j =1,2,...,8; /=l,2,...,ne

Копирование массива wa (k) = w(k) , k = 1,2,...,ng

Запись массивов eia (ne, 8) и eiw (ne, 8) на диск

Рис. 2. Блок-схема программы нелинейного решателя: ne , ng - соответственно числа элементов и степеней свободы

ансамбля КЭ; eia(ne,8), eiw(ne,8) - массивы интенсивности деформаций в точках интегрирования КЭ; w (ng),

wa (ng) - рабочие массивы для хранения узловых перемещений, принадлежащих смежным итерациям iter ;

max - максимальное число итераций (защита от зацикливания); S=10-8 - положительное малое число, служащее

критерием сходимости МППУ

Константы материала Е = 2,05-10 МПа, V = 0,29. Условные пределы пропорциональности и пределы временного сопротивления принимаем по справочным данным равными: для стали 40Х а пц = 400 МПа,

а в = 726 МПа; для стали А-Ш а пц = 314МПа, а в =

= 367 МПа. Заданным величинам а пц соответствуют

значения интенсивности деформаций е ^ = 0,0017

(сталь 40Х) и е ,т = 0,0015 (сталь А-Ш).

Z3

А

п9п

Рис. 3. Расчетная схема куба 0,01x0,01x0,01 м МПа

600 -500 400 300 200 100 0

о., МПа

250 200 150 100 50

I 1 1 1 1 I

I 1 1 1 1 I

0,005 0,015 0,025 е, а)

0 0,005

0,015

б)

0,025 е,

Рис. 4. Диаграммы деформирования о , ~ е , : а - сталь 40Х, п = 0,069; б - сталь А-Ш, п = 0,001

Аналитическое решение рассматриваемой статически определимой задачи имеет вид

E

f

1 -

1 -

е , ^

iT

(1 -п)

е = q , е > е

Значения интенсивности распределенной нагрузки принимались для куба из стали 40Х q = 500МПа, для куба из стали А-III - q = 365 МПа. Считаем, что нагрузка прикладывается мгновенно. Дискретизация куба в комплексах POLYGON и ABAQUS задавалась сетью 5x5x5 элементов в форме «кубиков». В комплексе COSMOS использовано автоматическое разбиение на 6852 элемента в виде тетраэдров. Анализ результатов сведен в табл. 1.

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

Пример 2. Расчет толстостенного цилиндра, нагруженного внутренним давлением q (задача о плоской деформации). Размеры сечения цилиндра: внутренний радиус r 1 = 33,5 мм, наружный радиус r 2 =

= 93,5 мм. Расчетная схема показана на рис. 5.

Z,

Рис. 5. Толстостенный цилиндр, загруженный внутренним давлением д

Фактически это тонкий слой (толщина 5 мм), вырезанный из цилиндра двумя параллельными поперечными сечениями, на которые наложены скользящие связи в осевом направлении. Константы материала (для стали) Е = 2-105 МПа, V = 0,29, е ,т = 0,002.

Схематизированная диаграмма деформирования представлена на рис. 6.

а,, МПа

500 400 300 т 200 -100 0

0,005

0,015 0,025 е,

Рис. 6. Диаграмма оt ~ е . для п = 0,034

Z

2

Таблица 1

Материал Решение МКЭ Аналитическое решение

POLYGON ABAQUS COSMOS

Сталь 40Х uZ3 = —0,1205 10-3 м uZl = 0,5511 •Ю-4 м uZ3 = -0,9464-10-4 м uZl = 0,4251 •Ю-4 м uZ3 = -0,1222-10-3 м uZl = 0,3629-10-4 м uZ3 = —0,127110-3 м

Сталь А-III uZ3 = -0,2495-10-2 м uZl = 0,1244 10-2 м uZ 3 = -0,221310-2 м uZl = 0,1328 10-2 м Решатель «зациклился» uZ3 = -0,2510-10-2 м

Сравнительные данные по перемещениям на внутреннем u r 1 и наружном u r радиусах цилиндра,

полученные в линейно упругой постановке (задача Ламе) для q = 500 МПа, приведены в табл. 2. В комплексах POLYGON и ABAQUS использовались объемные восьмиузловые КЭ соответственно на сетках 12x18 и 12x16 (1 слой элементов в осевом направлении), в комплексе COSMOS - автоматическое разбиение на 7892 элемента (тетраэдры). Установлено, что используемые схемы дискретизации дают результаты, достаточно хорошо совпадающие с точным решением.

Таблица 2

Комплекс Перемещение

url ur 2

POLYGON 0,1306-10-3м 0,6299^10-4м

ABAQUS 0,1306^10-Зм 0,5987^10-4м

COSMOS 0,1305^10-Зм 0,6294-10-4м

Аналитическое решение 0,1306^10-Зм 0,6306^10-4м

В проведенном числовом эксперименте выявлено, что уменьшением шага приращения нагрузки в арифметической прогрессии приводит к росту перемещений и г 1 и и г в упругопластической области работы

материала в геометрической прогрессии.

иг1-103, м 1,4-

H/м2

а)

ur240 , м 0,6

Данные сравнительного анализа с учетом пластического деформирования материала, для тестируемых коммерческих комплексов, приведены на рис. 7. Как видно комплекс COSMOS дает более близкие к точному решению результаты, чем комплекс ABAQUS. Вместе с тем расчетное время при использовании COSMOSa намного больше, чем ABAQUSa.

На рис. 7 также представлены графики u ~ q и

u r2 ~ q, полученные с помощью комплекса

POLYGON при различных способах управления процессом активного нагружения. Практический интерес представляет вариант аддитивного увеличения нагрузки (кривые 1, 2, 3), позволяющий опосредованно учитывать скорость деформирования. Общеизвестно, что повышение скорости силового воздействия приводит к понижению пластичности металла. Здесь в качестве шага по временной координате выступает величина приращения нагрузки Aq (чем больше величина Aq, тем выше скорость нагружения). Согласно полученным данным наращивание внутреннего давления в цилиндре от 0 до 500М Па за четыре шага ( Aq = 2,5-108 Н/м2) наиболее близко приближает ко-нечноэлементное решение к аналитическому (кривая 1).

H/м2

б)

Рис. 7. Графики перемещений на внутренней и и наружной un поверхностях цилиндра: 1 - qk = kAq , k = 1,2,...,4, (Aq =2,5-108Н/м2); 2-qk = kAq , k = 1,2,...,5 , (Aq =1-108 Н/м2); 3 -qk = kAq , k = 1,2,...,10 ,(Aq =0,5-108 Н/м2); 4 - Aqk =1-108Н/м2, k = 1,2,...,5 ; 5- Aqk = 0,5-108Н/м2, k = 1,2,...,10 ; 6- точное решение; 7 - решение COSMOS; 8 - решение ABAQUS

Исследования показали, что в случае дискретного нагружения цилиндра при изменении шага нагрузки Aq k в два раза перемещения u и u Г2 изменяются

незначительно (кривые 4, 5).

а)

V = 0,29, е ,т = 0,0015. Учитывая симметрию геометрии и нагрузки, в качестве расчетной схемы принимаем 1/4 часть пластины (рис. 9).

А

qq

Zl

б)

Рис. 8. Картины распределения радиальных напряжений о г в цилиндре при аддитивном нагружении = кДд ,

Дд = 2,5-108 Н/м2: а - к = 1,2; б - к = 1,2,...,4

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

Пример 3. Круглая пластина радиусом г = 0,2 м и толщиной к = 0,01 м, жестко закрепленная по контуру и нагруженная односторонним давлением д . Материал пластины - сталь с константами Е = 2,05-105 МПа,

Z,

Рис. 9. Схема разбивки на КЭ (комплекс POLYGON) и закрепления 1/4 части круглой пластины

Для оценки качества дискретизации пластины объемными КЭ сначала выполнялся расчет в линейно упругой постановке для q =1 МПа. Данные теста

представлены в табл. 3. Следует отметить, что в «банке» комплекса ABAQUS предусмотрены специальные пластинчатые элементы типа shell. Однако для установления равных стартовых условий для тестируемых программ специально выбраны объемные КЭ типа solid и активизирована опция «большие перемещения» (576 восьмиузловых КЭ). В исследовательском комплексе POLYGON использованы объемные КЭ двух типов - восьмиузловые и шестиузловые, примыкающие к полюсу пластины (всего 512 КЭ). Расчеты в POLYGONе и ABAQUSe выполнялись с применением двух слоев КЭ по толщине пластины, что улучшает конечно-элементное решение при v ф 0 по сравнению с однослойной схемой разбивки. В комплексе COSMOS генерация конечно-элементной сетки, как и в предыдущих примерах, осуществлялась автоматически (задействовано 1107 КЭ в виде тетраэдров).

Таблица 3

Комплекс Перемещение в центре пластины UZ 3max ' м

POLYGON 0Д247-10-2

ABAQUS 0Д780-10-2

COSMOS 0,1336 10-2

Аналитическое решение 0,1340 10-2

Результаты расчетов пластины в виде значений перемещений в центре пластины и х тах и графиков

Z 3 max

рис. 10.

- q с учетом пластичности показаны на

Uz3 max, м

4 q-10-6, H/м2

Рис. 10. Графики перемещений в центре круглой пластины: 1- qk = kAq , к = 1,2,...,20 , ( Aq =0,25-106 Н/м2); 2 - qk = kAq , к = 1,2,...,10, ( Aq =0,5-106 Н/м2); 3 - qk = kAq , k = 1,2,...,5 , ( Aq = 1,0106 Н/м2); 4 - qk = kAq , k = 1,2, ( Aq =2,0-106 Н/м2);

5 - решение ABAQUS и COSMOS (q = 3,5-106 Н/м2); 6 - решение POLYGON (q = 3,5-106 Н/м2)

Из представленных данных видно, что при мгновенном приложении давления q = 3,5-106 Н/м2 решения в COSMOSе и ABAQUSе совпадают. Программа POLYGON в этом случае дает заниженное значение перемещения uZ 3max. Далее исследовалось влияние величины шага приращения давления A q на прогиб uZ max при аддитивном законе

нагружения. Установлено, что при A q = 2-106 Н/м2 (кривая 4) значение u Z 3 max , полученное с помощью

POLYGONa, приближается к результатам COSMOSa и ABAQUSa.

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

0.00080541 0.00161082 0.00241 Б2Э

0.00322164

б)

Рис. 11. Визуализация деформации круглой пластины (масштаб 1:1): а - 9 = 2,5-106 Н/м2; б - 9 = 5,0106 Н/м2

Литература

1. Гайджуров П.П. Конечно-элементное решение упруго-пластической задачи при циклическом нагружении // Изв. вузов. Сев.-Кавк. регион. Технические науки. -2003. - Приложение № 3.- С. 83-87.

2. Биргер И.А., Мавлютов Р.Р. Сопротивление материалов. -М.: Наука, 1986.-560 с.

3. Лукаш П.А. Основы нелинейной строительной механики. -М.: Стройиздат, 1978. - 204 с.

4. Ильюшин А.А. Пластичность. - М.: ОГИЗ, ГИТТЛ, 1948.-

376 с.

Южно-Российский государственный технический университет (Новочеркасский политехнический институт)

11 сентября 2006 г.