CC BY
34
Известия Тульского государственного университета Естественные науки. 2013. Вып. 1. С. 99-110
Механика =
УДК 539.3
Конечноэлементный анализ предельного формоизменения тонкого алюминиевого листа при осесимметричном гидровыпучивании в матрицу с плоским дном
В. Л. Михайлова, В. К. Петров, Л. Г. Сухомлинов
Аннотация. Излагаются результаты применения осесимметричной жесткопластической безмоментной конечноэлементной модели к исследованию предельного формоизменения тонкого алюминиевого листа при гидровыпучивании в матрицу с плоским дном. В процессе численного моделирования ситуация разрыва формуемой оболочки определяется на основе сравнения деформаций в каждом из элементов модели с диаграммой предельных деформаций исследуемого листового материала.
Ключевые слова: осесимметричная жесткопластическая безмо-ментная конечноэлементная модель, алюминиевый лист, гидровыпучивание в матрицу, предельное формоизменение.
К настоящему времени опубликовано множество результатов расчетноэкспериментальных исследований в подтверждение способности известных конечноэлементных комплексов, ориентированных на исследование процессов пластического формоизменения листовых металлов, давать достоверный прогноз по силовым и деформационным характеристикам подобных процессов. В гораздо меньшей степени представлены расчетно-эксперименальные результаты, демонстрирующие способность тех же вычислительных комплексов давать надежный прогноз по моменту разрыва формуемого листового металла в исследуемом процессе формоизменения. Подобная ситуация во многом обусловлена невысокой степенью точности используемых при таких прогнозах диаграмм предельных деформаций листовых металлов, определяемых на основе экспериментов. В последнее время, в связи с совершенствованием техники эксперимента, стали появляться уточненные данные по упомянутым диаграммам для ряда листовых металлов. В частности, в работах [1-3], приведены подобные уточненные данные, касающиеся тонкого листового алюминия (чистый алюминий: А199,5; AW1050). Излагаемое ниже
имеет целью на примере процесса формовки тонкого алюминиевого листа при осесимметричном гидровыпучивании в матрицу с плоским дном (экспериментально исследованном в работе [4]) продемонстрировать возможности осесимметричной жесткопластической безмоментной конечноэлементной модели [5,6] по предсказанию поведения формуемой оболочки и определению момента ее разрыва (с учетом результатов работ [1-3]).
Основные положения используемой при исследованиях вычислительной модели состоят в следующем.
Исходим из предположения, что формуемая из листового металла под действием давления рабочей жидкости и жестких инструментов осесимметричная оболочка относится к классу тонких безмоментных оболочек. Упругими деформациями на фоне больших пластических деформаций пренебрегаем, считая материал оболочки жесткопластическим. Используем предложенный Р. Хиллом [7] вариант теории течения (квадратичный критерий текучести) для трансверсально изотропного материала с изотропным упрочнением. Считаем, что взаимодействие оболочки с инструментом осуществляется в соответствии с кулоновским законом трения. Меридиан срединной поверхности рассматриваемой оболочки в ее исходном недеформированном состоянии разбиваем на такое количество N участков малых размеров, чтобы в течение всего процесса деформирования допустимо было бы пренебрегать их кривизной, считая эти участки прямолинейными. Процесс формоизменения подобной безмоментной оболочечной модели, состоящей из указанных N элементарных оболочек с прямолинейными образующими, рассматриваем как пошаговый, при котором переход из известного состояния в момент времени Ь в новое состояние, относящееся к моменту времени Ь + АЬ, осуществляется с малыми приращениями деформаций. На указанном малом временном интервале АЬ (шаге нагружения) формулировку задачи для принятой дискретной модели оболочки выполняем в терминах узловых перемещений (перемещений концов отрезков разбиения упомянутого меридиана оболочки) с использованием цилиндрической системы координат (х, г, ф).
Обозначая через х, г и х*, г* координаты точек срединной поверхности дискретной модели рассматриваемой оболочки в начале и конце текущего шага нагружения, а через их, иг осевые и радиальные перемещения указанных точек на данном шаге, мы можем записать
х* = х + их, (х ^ г). (1)
Координаты и перемещения, относящиеся к узлам дискретной модели, будем обозначать как хг, гг, х*, г*, игх, игг (г = 1, 2,N + 1). С учетом этого, принимая во внимание связи (1), для значений в и в* длины прямолинейной образующей г-го элемента (г = 1, 2,N) в начале и конце шага нагружения, можно записать
в = V(хг+1 - хг)2 + (Гг+1 - Гг)2,
в* = \! (х*+1- х*)2+(г*+1- г*)2 = (2)
= \/(хг+1 - хг + иХ+1 - иХ)2 + (гг+1 - гг + иГ+1 - игг)2.
Здесь и далее индекс г у величин, относящихся к г-му элементу, для простоты опускается.
Приращение деформации (малое относительное удлинение) элемента модели на шаге нагружения в меридиональном направлении (е8) определяем с учетом (2) по схеме
в * _ в в *2 — в2 1
ез = ---в-- ~ ---2в---- = 2-2 [2(х*+1 - хг)(и%х+1 - иХ) +
в 2в 2в (3)
+2(гг+1 - гг)(игг+1 - игг) + (игх+1 - игх)2 + (игг+1 - игг)2 ].
Отмечаем, что выражение (3) является квадратично нелинейным по отношению к узловым перемещениям на шаге нагружения. Оно справедливо для любых значений этих перемещений при условии, что |е8| ^ 1.
Малое относительное удлинение на шаге нагружения окружного материального волокна в середине элемента (приращение деформации еф в окружном направлении в середине элемента) определяем с учетом (1) по схеме
еФ = (г* - г)/г = иг/г, (4)
где
г = (гг + гг+1)/2, Г * = (г* + г*+1)/2, и = (игг + игг+1 )/2. (5)
Приращение деформации еь по толщине элемента на шаге нагружения
определяем, исходя из условия несжимаемости, а именно
еН = -(ев + еф)- (6)
Основываясь на принятом варианте теории течения, соотношения между устанавливающимися в элементах модели на шаге нагружения напряжениями и соответствующими приращениями деформаций записываем в виде
Ае = (1+ 2Я) ^(еФ - ^ен)2 + (Яен - ев)2 + Я(ея - е^2 (7)
а = Ф(е + Ае),
где аэ и аф — меридиональное и окружное напряжения; а — эквивалентное напряжение; Ае — приращение накопленной эквивалентной деформации; е и е + Ае = е* — величины накопленной эквивалентной деформации в начале и конце шага нагружения; Я и а = Ф(е) — экспериментально определяемые параметр нормальной анизотропии и функция упрочнения листового металла.
Вводим обозначение д для величины давления, приложенного к поверхности каждого элемента модели и действующего в направлении нормали к этой поверхности, а также обозначения 2пГХ и 2п¥гг (г = 1, 2, ...,Ы + 1) для обобщенных сил в радиальном и окружном направлениях, приложенных к узлам модели. Исходим из принципа возможных перемещений, утверждающего, что в состоянии равновесия работа приложенных к деформируемому телу сил на вариациях возможных перемещений равна работе напряжений на соответствующих вариациях деформаций. В результате для рассматриваемой модели оболочки на шаге нагружения получаем вариационное уравнение вида
2п^2 (аэ $ез + аф беф) гкв = 2п^д (п*х бь,х + п* би) гв+
М+1 (8)
+2п ^2 (Гх;бигх + Еггбигг). (
г=1
Суммирования без указания индексов в уравнении (8) осуществляются по всем элементам модели. Величины пх и п* представляют собой осевую и радиальную компоненты единичного вектора нормали к поверхности элемента, определяемые по схеме
п* = (х*+1 - х*)/в* & (хг+1 - хг + игх+1 - игх)/-,
(9)
п*х = -(г*+1 - г*)/в * & -(гг+1 - гг + К+1 - игг)/в.
С использованием связей (3)-(7), (9) вариационное уравнение (8) формулируется исключительно в терминах узловых перемещений. Приравнивая коэффициенты при одноименных вариациях узловых перемещений в левой и правой части такого уравнения, получаем для г-го узла дискретной модели оболочки на шаге нагружения уравнения равновесия вида
ьх(их,и3г) = РхЦ, Ьгг(их,и3г) = ¥1, (] = г - 1, г, г + 1), (10)
где Ьгх, Ьгг — нелинейные алгебраические операторы, отражающие геометрическую и физическую нелинейность рассматриваемой задачи на шаге нагружения.
В случае закрепленного г-го узла относящаяся к этому узлу пара силовых уравнений (10) заменяется кинематическими условиями вида
игх = 0, игг = 0.
(11)
Если г-й узел находится в свободном состоянии (не контактирует с инструментом), в соответствующих ему уравнениях (10) следует положить
¥х = 0, Рг = 0. (12)
Для узла, находящегося в контакте с инструментом, контур которого на шаге нагружения описывается уравнением ^(х,г) = 0, должно выполняться кинематическое условие вида
^(хг + игх,гг + игг ) = 0 (13)
и силовое условие (соответствующее кулоновскому закону трения) вида
I -1 Аиг
К = -Ц ¥I , (14)
е р 1 п \Аиге\ ’ у 7
где ц — коэффициент трения; Аиге — относительное перемещение узла вдоль поверхности инструмента на шаге нагружения; ¥1 и — тангенциальная и
нормальная по отношению к контуру инструмента компоненты вектора узловой силы. Обозначая как ех, егг и пгх, пгг компоненты векторов касательной и нормали к контуру инструмента в текущем положении г-го узла, силы К, , входящие в уравнение (14), можно выразить через силовые факторы ¥х, ¥£ в виде
К = екх + ег , ¥гп = пгх¥х + пГ ¥гг. (15)
Использование записи (15) предполагает, что входящие в нее факторы
¥х, посредством связей (10) явно выражены через искомые узловые пере-
мещения на шаге нагружения.
Что касается перемещения Аиге, входящего в формулировку уравнения (14), то для него имеет место зависимость вида
Аиге = егх(игх - Аи) + еггигг, (16)
где Аи — заданное перемещение инструмента в направлении оси х на шаге нагружения.
Замечаем, что физические соотношения (7) и силовые соотношения (14) имеют особенности при Ае = 0 и 1А^^1е1 = 0. В целях устранения этих особенностей предусмотрена регуляризация указанных соотношений путем введения в пренебрежимо малых окрестностях значений Ае = 0 и |Аие| = 0 соответствующих фиктивных линейных относительно Ае и | Аие| участков. Подобный прием [5, 6] позволяет напряжениям аэ, аф и силам ¥'I при необходимости снижаться вплоть до нулевых значений, обеспечивая выполнение условий равновесия на шаге нагружения.
Решение физически и геометрически нелинейной контактной задачи для описанной дискретной модели оболочки на шаге нагружения сводится посредством итерационной процедуры к решению последовательности линей-
ных задач. При этом линеаризация исходной нелинейной системы уравнений на шаге нагружения в рамках такой процедуры осуществляется с использованием методов Ньютона и переменных параметров. Итерационные уточнения продолжаются до достижения заданной относительной точности (бот) по перемещениям. Решение соответствующей системы линейных алгебраических уравнений выполняется по методу Гаусса.
Перейдем теперь к заявленному исследованию с применением описанной модели.
На рис. 1 представлена схема гидровыпучивания (под действием давления д) круговой заготовки из тонкого алюминиевого листа (толщиной Н = 0, 31мм) в цилиндрическую матрицу (радиусом а = 50мм) с плоским дном. Переферийная зона заготовки жестко закреплена. С учетом малости радиуса скругления рабочей кромки матрицы (составляющего величину порядка 2мм) моделирование осуществляем применительно к случаю круговой заготовки, закрепленной по контуру радиуса а = 50мм.
Согласно данным работы [4], рассматриваемый листовой алюминий представляет собой изотропный материал (Я = 1). Из испытаний на одноосное растяжение была установлена кривая упрочнения данного материала в виде а = Леа (где А = 156, 4МПа, п = 0, 29). Значение коэффициента трения (ц = 0, 3) в зоне контакта формуемого листа с пластиной, представляющей собой дно матрицы, было оценено в испытаниях с протягиванием по поверхности указанной пластины образцов из исследуемого листового алюминия, прижимаемых к пластине заданной вертикальной нагрузкой. Эксперименты по гидровыпучиванию описанной листовой заготовки проводились с использованием матриц глубиной Ь = 10мм и Ь = 15мм. В процессе гидровыпучивания давление д постепенно увеличивалось от нуля до 1,4МПа. При этом в случае Ь = 10мм разрыва формуемой оболочки не наблюдалось, а в случае Ь = 15мм такой разрыв был зафиксирован при достижении давлением д значения 1,4МПа.
При численном моделировании в качестве параметра нагружения принималось давление д, изменяемое с шагом Ад = 0, 001МПа. Кроме того было принято N = 200 и бот = 0, 001.
На рис. 2 представлены относящиеся к случаю Ь = 10мм расчетные (сплошные кривые) и экспериментальные (точки) результаты по распределению меридиональной логарифмической деформации еэ вдоль радиуса г рассматриваемой круговой заготовки при д = 1, 38МПа. Кривые, помеченные цифрами 1, 2, 3, 4, 5, 6, 7, соответствуют вариантам расчетов с выбором значений параметров ц, А [МПа], п в виде (0,3; 156,4; 0,29), (0,4; 170,0; 0,25), (0,4; 156,4; 0,25), (0,4; 156,4; 0,29), (0,43; 156,4; 0,29), (0,47; 156,4; 0,29), (0,4; 156,4; 0,32). Видно, что задание коэффициента трения в виде ц = 0, 3 (в соответствии с оценкой, данной в работе [4]) приводит к заметно завышенным по сравнению с экспериментом значениям деформаций исследуемой оболочки в зоне ее контакта с дном матрицы. К более реалистичной картине в этом плане приводят значения ц порядка 0,4. Отмечаем также, что изменение
Рис. 1. Схема гидровыпучивания круговой заготовки из тонкого алюминиевого листа в матрицу с плоским дном
значений параметров А и п кривой упрочнения рассматриваемого листового алюминия в ту или иную сторону относительно значений А = 156, 4МПа и п = 0, 29 (установленных в испытаниях на одноосное растяжение) ведет к заметному отклонению расчетных результатов от эксперимента. Исходя из этого, при проведении расчетных исследований по прогнозу момента разрыва формуемой оболочки принимаем л = 0, 43, А = 156, 4МПа, п = 0, 29.
4-С> г^:5/6
%/ /
1-" Л 2
Рис. 2. Графики распределения меридиональной логарифмической деформации вдоль радиуса рассматриваемой круговой заготовки при различных значениях параметров л, А, п в случае Ь = 10мм и
ц =1, 38МПа
На рис. 3, 4 представлены относящиеся к случаю матрицы глубиной Ь = 15мм результаты по распределению логарифмических меридиональной е3 и окружной еф деформаций вдоль радиуса г рассматриваемой круговой
заготовки при различных значениях давления ц. Расчетные кривые, помеченные цифрами 1,2,3,4,5, соответствуют значениям давления ц = 1,26; 1,32; 1,4; 1,5; 1,6МПа. Точками изображены экспериментальные результаты, относящиеся к ситуации ц = 1, 32МПа. Отмечаем хорошее согласование расчетных и экспериментальных результатов по распределению деформаций в данной ситуации. Замечаем, что в исследуемом процессе формоизменения при значениях давления 1,4МПа и выше (варианты 3,4,5) у соответствующих кривых распределения деформаций наблюдается характерный всплеск, охватывающий небольшой по протяженности участок формуемой оболочки вблизи зоны контакта с дном матрицы. С ростом давления ц уровень деформаций на указанном участке существенно возрастает, грозя привести к разрыву формуемого листа.
Здесь следует обратить внимание на то обстоятельство, что разрыв таких высокопластичных материалов, каким является рассматриваемый листовой алюминий, связан с таким явлением, как локализация деформации. Суть этого явления состоит в том, что на каком-то этапе процесса пластического формоизменения в каком-то из мест растягиваемого участка формуемого листового материала (под влиянием всегда присутствующих в материале микронеоднородностей) образуется локальное утонение («шейка»). При этом рост деформации в зоне, примыкающей к области «шейки», прекращается. А продолжающее (катастрофически) развиваться локальное утонение приводит к разрыву листового материала вдоль линии, проходящей через область «шейки». С учетом этого, предельные значения деформаций (значения деформаций в момент образования «шейки») для подобных листовых материалов определяют из испытаний соответствующих образцов, растягивая их до разрыва и замеряя деформации в примыкающих к линии разрыва элементах таких образцов. Вместо описанных экспериментов с образцами используют и вычислительные подходы, при которых решается задача теории пластичности для растягиваемого тонкого листа, имеющего пренебрежимо малую начальную неоднородность (например, в форме узкого участка с незначительно уменьшенной начальной толщиной). Предельные значения деформаций определяют при этом, фиксируя момент начала локального развития деформации.
На рис. 5 треугольниками отмечены данные (взятые из статей [1, 2]) по предельным деформациям листового алюминия, полученные из испытаний алюминиевых образцов толщиной Н = 1мм (темные треугольники) и Н = 0,1мм (светлые треугольники). Здесь и £2 — предшествующие разрыву значения главных логарифмических деформаций (максимальной и минимальной) в элементах исследуемого листового материала вблизи линии разрыва. Как видно, при £2 > 0 с увеличением толщины алюминиевого листа от 0,1 до 1мм предельно допустимые для его безразрывного формоизменения значения деформаций несколько возрастают. На такую тенденцию указывают также и результаты экспериментальных исследований, представленные в работе [3].
А- 5
1:
-• *" • 1
о » 20 зо 40 г, мы
Рис. 3. Графики распределения меридиональной логарифмической деформации вдоль радиуса рассматриваемой круговой заготовки при различных значениях давления д в случае Ь = 15мм
/5
\
2 ' %
Л
\
•\
О 10 20 30 40 }', ЛШ
Рис. 4. Графики распределения окружной логарифмической деформации вдоль радиуса рассматриваемой круговой заготовки при различных значениях давления д в случае Ь = 15мм
% ч % *. * ^
0.2 -0.1 0 0.1 0.2 8^
Рис. 5. Диаграммы предельных деформаций образцов листового алюминия с толщинами Н = 1мм (темные треугольники) и Н = 0,1мм (светлые треугольники)
Это дает основание полагать, что диаграмма предельных деформаций для интересующего нас случая листового алюминия толщиной h = 0, 31мм занимает при £ 2 > 0 некоторое промежуточное положение между диаграммами, относящимися к случаям h = 0,1 и h = 1мм. На рис. 5 приведена также (пунктиром) кривая предельных деформаций для листового алюминия, полученная в работе [1] указанным выше расчетным путем. Учитывая, что при £2 > 0 она занимает упомянутое промежуточное положение между двумя экспериментально установленными диаграммами, эту кривую и примем в качестве диаграммы предельных деформаций для исследуемого листового алюминия.
Снова обращаемся к результатам, представленным на рис. 3, 4, относящимся к случаю матрицы глубиной b = 15мм. Определяя уровень деформаций £s и £ф в наиболее деформированной зоне формуемой оболочки, отмечаем, что при q = 1, 6МПа это величины порядка 0,2 и 0,11 , а при q = 1, 5МПа — порядка 0,17 и 0,11. Замечая, что точка с координатами £\ = 0, 2 и £2 = 0,11 (см. рис. 5) расположена выше принятой кривой предельных деформаций, а точка с координатами £\ = 0,17 и £2 = 0,11 практически попадает на эту кривую, приходим к выводу, что разрыв формуемой оболочки следует ожидать при давлении q порядка 1,5МПа. Указанный прогноз хорошо согласуется с экспериментально установленным фактом разрыва оболочки при q = 1, 4МПа.
В качестве общего вывода по изложенной статье отметим, что представленные (и подтвержденные сравнением с экспериментом) результаты расчетных исследований позволили на примере процесса гидровыпучивания тонкого алюминиевого листа в матрицу с плоским дном продемонстрировать возможности применения осесимметричной жесткопластической безмо-ментной конечноэлементной модели как в предсказании деформационных характеристик, так и в прогнозе моментов разрыва формуемых оболочек в процессах формоизменения листовых металлов.
Список литературы
1. Ogorek A., Stachowicz F. Determination of forming limits of thin aluminium sheets // Int. Multidisciplinary Conf. 6-th edition / Baia Mare. 2005. P. 551-556.
2. Diehl A., Vierzigmann U., Engel U. Characterisation of the mechanical behavior and the forming limits of metal foils using a pneumatic bulge test // Int. J. Mater. Form. 2009. V. 2. N 1. P. 605-608.
3. Diehl A., Staud D., Engel U. Investigation of the mechanical behaviour of thin metal sheets using the hydraulic bulge test // Proceedings of the 4th Int. Conf. on Multi-Material Micro Manufacture / Dunbeath, Scotland, UK: Whittles Publishing, 2008. P. 195-198.
4. Nakamachi E., Takezono S., Sowerby R. A numerical analysis of the hydraulic bulging of circular disks into axisymmetric dies // Trans. ASME. J. Appl. Mech. 1982. V. 49. N 3. P. 501-506.
5. Сухомлинов Л.Г., Энгельсберг В.К. Конечноэлементная система автоматизированного расчета напряженно-деформированного состояния тонких оболочек в процессах осесимметричного формоизменения под действием жестких штампов // Изв. вузов. Машиностроение. 1989. N 3. С. 66-71.
6. Sukhomlinov L.G., Engelsberg V.K., Davydov V.N. A finite element membrane model for the analysis of axisymmetric sheet metal forming processes // Int. J. Mech. Sci. 1992. V. 34. N 3. P. 179-193.
7. Хилл Р. Математическая теория пластичности. М.: ГИТТЛ, 1956. 407 с.
Михайлова Виктория Львовна (tm@mami.ru), к.т.н., доцент, кафедра теоретической механики, Московский государственный машиностроительный университет (МАМИ).
Петров Владимир Кириллович (tm@mami.ru), к.т.н., доцент, кафедра теоретической механики, Московский государственный машиностроительный университет (МАМИ).
Сухомлинов Лев Георгиевич (tm@mami.ru), д.т.н., профессор, кафедра теоретической механики, Московский государственный машиностроительный университет (МАМИ).
A finite element analysis of forming limit of a thin aluminium
sheet being bulged by hydraulic pressure into an axisymmetric
flat bottomed die
V. L. Mikhaylova, V. K. Petrov, L. G. Sukhomlinov
Abstract. The results are presented on the application of an axisymmetric rigid-plastic finite element membrane model to the analysis of forming limit of a thin aluminium sheet being bulged by hydraulic pressure into a flat bottomed die. In the process of numerical simulations, the situation with failure of the deformed shell is recorded on the basis of comparing the computed strains in each membrane element with the forming limit diagram of the sheet material under investigation.
Keywords: axisymmetric rigid-plastic finite element membrane model, aluminium sheet, hydraulic bulging into die, forming limit.
Mikhaylova VKtoriya (tm@mami.ru), candidate of technical sciences, associated professor, department of theoretical mechanics, Moscow State University of Mechanical Engineering (UMech).
Petrov Vladimir (tm@mami.ru), candidate of technical sciences, associated professor, department of theoretical mechanics, Moscow State University of Mechanical Engineering (UMech).
Sukhomlinov Lev (tm@mami.ru), doctor of technical sciences, professor, department of theoretical mechanics, Moscow State University of Mechanical Engineering (UMech).
Поступила 11.12.2012
