Математична модель процесу витяжки трубчастої заготовки
1. Варіаційні підходи до вирішення задач методом кінцевого елемента
Основна ідея МКЕ грунтується на заміні деякої безперервної величини в межах даної області дискретною моделлю, яка будується на безлічі кусково-неперервних функцій, визначених на кінцевому числі підобластей, званих кінцевими елементами (КЕ). Невідома шукана величина в межах кожного КЕ апроксимується, як правило, поліноміальної функцією заданого виду з урахуванням вимоги безперервності на кордонах суміжних КЕ. При цьому вибір форми кінцевого елементу та виду вираження, що апроксимує дійсний закон зміни досліджуваної величини в межах КЕ, є одним з найбільш відповідальних моментів у загальній процедурі МСЕ, від якого істотно залежить точність наближеного рішення. Таким чином, безперервна в межах досліджуваної області невідома величина (наприклад, переміщення, швидкість переміщення, напруга, температура і т. д.) представляється через кінцеве число її дискретних значень у вузлах елементів.
Побудова дозволяють рівнянь МКЕ для вирішення задач механіки деформівних середовищ базується на відповідних варіаційних принципах і випливає з оптимізації деякої інтегральної величини (функціоналу), пов'язаної з роботою або потужністю напруг і зовнішньої прикладеного навантаження при дотриманні заданих граничних умов. У загальному вигляді такий функціонал з урахуванням дії масових і поверхневих сил можна представити виразом:
, (7)
де A Д - робота або потужність внутрішніх сил; A М - робота або потужність, що розвивається масовими силами; A В - робота або потужність зовнішніх сил.
Подальша процедура МКЕ передбачає подання вираження (2.1) у вигляді функціонала значень, невідомих тільки у вузлах КЕ, і побудова роздільної системи рівнянь шляхом мінімізації J по всіх вузловим змінним:
(7)
Проте, зазначений спосіб отримання дозвільних рівнянь для КЕ за допомогою функціонала (1) не є єдино можливим. В даний час рівняння для елементів отримують шляхом мінімізації функціонала, пов'язаного з даним диференціальним рівнянням відповідної задачі математичної фізики. Відомі також кінцево-елементні рішення, засновані на методі Гальоркіна. В останньому випадку відпадає необхідність у варіаційної формулюванні задачі.
Спосіб отримання дозвільних рівнянь для КЕ, заснований на оптимізації функціоналу (1), є загальновизнаним при теоретичному вирішенні завдань ОМД, оскільки варіаційні принципи мають наочний фізичний зміст і досить суворе математичне обгрунтування.
По відношенню до функціонала (1) відомі три види варіаційних принципів теорії пластичності в залежності від того, через які змінні величини виражена потужність (потенційна енергія) деформації.
Принцип мінімуму повної потужності (повної енергії) чи принцип можливих змін деформованого стану розглядає потужність (потенційну енергію) деформівного тіла як функціонал довільної системи швидкостей (переміщень), що задовольняє кинематическим граничним умовам, і який приймає мінімальне значення для системи швидкостей (переміщень), фактично реалізованої в деформується тілі.
Принцип мінімуму додаткової роботи Кастільяна або принцип можливих змін напруженого стану розглядає додаткову роботу як функціонал довільної системи напруг, що задовольняє рівнянь рівноваги усередині тіла і на його поверхні, і який приймає мінімальне значення для системи напруг, фактично реалізованої в деформується тілі.
У варіаційному принципі Рейсснера або принципі можливих змін напруженого і деформованого станів потужність (енергія) розглядається як функціонал швидкостей і напруг, і змінні тієї й іншої групи варіюються незалежно один від одного.
Кожному з перерахованих варіаційних принципів відповідає певна форма МСЕ. Принципу мінімуму повної потужності (повної енергії) відповідає кінематичний метод, принципом мінімуму додаткової роботи - метод напружень, а варіаційного принципом Рейсснера - змішаний метод.
При навантаженні тіла потенційна енергія зовнішніх сил змінюється. При цьому зовнішні сили роблять роботу. Потенціал зовнішніх сил Q на можливих переміщеннях δu чисельно дорівнює роботі цих сил:
(7)
де P i - поверхневі сили, S - площа поверхні тіла.
У результаті зміни потенційної енергії зовнішніх сил тіло деформується і накопичує потенційну енергію деформації W
(7)
де s ij - компоненти тензора напруги, e ij - компоненти тензора деформації, V - об'єм тіла.
Сума енергії деформації і потенціалу зовнішніх сил дорівнює повній потенційної енергії:
(7)
Відповідно до принципу можливих переміщень Лагранжа зміна повної потенційної енергії на можливих переміщеннях дорівнює нулю:
(7)
При цьому під можливими переміщеннями d u розуміються як завгодно малі відхилення системи від положення рівноваги, що допускаються накладеними на систему зв'язками. З рівняння (6) випливає, що в стані рівноваги енергія П має стаціонарне значення. Можна показати, що в положенні стійкої рівноваги цей екстремум відповідає мінімуму.
З урахуванням викладеного варіаційне рівняння Лагранжа для статичної задачі має вигляд:
(7)
Мінімізуючи потенційну енергію за можливим переміщенням, отримуємо систему лінійних рівнянь, розв'язуючи яку визначаємо значення зовнішніх сил.
2. Основні співвідношення методу скінченних елементів
Найпростішим елементом, застосовуваним для вирішення осесиметричної задачі механіки деформівного твердого тіла, є тороїдальний елемент із трьома вузлами, розташованими у вершинах трикутного перерізу (Рис. 1.).
Рис. 1. Кінцевий елемент в задачі осесиметричної деформації
Вектор переміщень вузлових точок кінцевого елемента у випадку осесиметричної деформації має вигляд:
. (8)
Довільна точка елемента отримує переміщення u r і u z у напрямку осей r і z. Тому матриця u має вигляд:
. (9)
Вузлові переміщення і u пов'язані між собою матрицею апроксимуючих функцій N:
(9 ')
Найбільш поширений спосіб отримання наближених рішень на основі використання варіаційного рівняння за методом Релея - Рітца. Він полягає в тому, що функції переміщень задаються у вигляді інтерполяційного полінома. Якщо обмежитися поліномом першого ступеня, то ці функції будуть мати вигляд:
(10)
Тут a i - довільні постійні. При лінійної апроксимації сторони трикутника після деформування елемента залишаються прямими.
Висловимо a i через переміщення вузлів елемента. У результаті матриця N прийме вигляд:
(11)
S - площа перерізу елемента:
, (12)
де r i, z i - координати i-го вузла у відповідних осях.
Деформований стан у будь-якій частині тіла описується тензором малих деформації Коші:
(13)
В умовах осесиметричної задачі тензор деформації другого рангу зводиться до вектора:
(14)
компоненти якої виражаються через похідні переміщень за відповідними координатами:
. (15)
Зв'язок між складовими векторів деформації і переміщень можна представити одним матричним рівністю:
(16)
де B - матричний диференційний оператор:
. (17)
Використовуючи (16) і (17), можна висловити деформації через вузлові переміщення
. (18)
Матриця функцій форми C для осесиметричної деформації:
. (19)
Коефіцієнти матриці C залежать від координат r і z точки всередині елемента. Для трикутника з клунками в вершинах координати r і z можна замінити середніми по елементу значеннями:
(20)
Вектор напружень s має вигляд:
(21)
Висловимо за допомогою лінійного закону, що виражається матрицею жорсткості, напруги через вузлові переміщення
, (21 ')
де D - матриця матеріальних констант.
Потенційна енергія деформації елемента з урахуванням (20) і (19)
. (22)
Інтеграл у вираженні (2.22) є матриця жорсткості вибраного елемента
, (23)
Елементарний обсяг . Тому матриця жорсткості елемента записується таким чином:
, (24)
де S - площа елемента.
З урахуванням виконаних перетворень рівняння рівноваги елемента через вузлові переміщення виражається у формі:
(25)
де K - матриця жорсткості; P, - Вектори зовнішніх сил і вузлових переміщень, відповідно.
При наявності пружних і пластичних деформації зв'язок між напруженнями і деформаціями нелінійна. Рішення нелінійної системи рівнянь дуже занадто. Тому при використанні деформаційної теорії часто використовують кусково-лінійний закон зв'язку напруг і деформації. Тоді при вирішенні завдання в збільшеннях напруг Ds і деформації De, зв'язок між якими можна вважати лінійної, отримуємо систему лінійних рівнянь:
(26)
Одним із способів вирішення завдання в збільшеннях є метод послідовних навантажень. Для квазистатическом завдання збільшення зовнішніх сил D P обчислюються на кроці за часом D t. При цьому вектор зовнішніх сил P в момент часу t дорівнює:
(27)
де n - крок навантаження.
Таким чином, з урахуванням вищевикладеного, варіаційне рівняння рівноваги в матричній запису приймає вигляд:
(28)
де - Вектор збільшень переміщень.
3. Представлення матриці жорсткості
У межах пружності зв'язок між приростами напруг і деформації виражається законом Гука. Згідно з ним компоненти збільшень деформації є лінійними функціями збільшень напруг. Пластичне стан матеріалу описується теорією малих пружно деформації Іллюшина. Приймається теорія ізотропного зміцнення. Об'ємна деформація у пластичній зоні залишається пружною і для неї виконується об'ємний закон Гука:
, (29)
де q - відносна зміна обсягу.
Модуль об'ємного стиснення k для ізотропного тіла у випадку осесиметричної деформації має вигляд:
. (30)
Модуль зсуву G пов'язаний з модулем Юнга E і коефіцієнтом Пуассона n формулою:
в пружній області:
(31)
у пластичній:
(32)
Тут H - дотичний модуль зміцнення. Коефіцієнт Ляме - l визначається формулою:
(33)
Таким чином, матриця матеріальних констант D має вигляд:
. (34)
Слід особливо відзначити, що використовувати матрицю жорсткості в такому вигляді для пластичного стану можна, тільки пов'язуючи збільшення деформації та напружень, про що було сказано раніше при виводі рівняння рівноваги.
Знаючи поточний стан елемента, межа плинності, накопичену деформацію і збільшення зовнішніх сил, можна визначити зміну напружено-деформованого стану на кроці збільшення переміщень D u і сил D Р, використовуючи для обчислення K за формулою пружне або пластичне подання матриці жорсткості.
4. Пластична деформація
Пластична деформація твердого тіла розглядається в рамках деформаційної теорії пластичності. Прийнято наступні вихідні положення:
тіло изотропно;
відносна зміна обсягу мало і є пружною деформацією, пропорційній середньому тиску:
або
;
повні збільшення складових деформації De ij складаються з збільшень складових пружною деформації De eij і пластичної деформації De pij:
;
девіатора збільшень напруги та деформації пропорційні:
.
Напружено-деформований стан елемента на i +1 кроці характеризується інтенсивністю деформації e i:
(35)
де e ij - компоненти тензора деформації.
Якщо інтенсивність деформації будь - якого кінцевого елемента перевищила поточний межа пружності по деформаціям , То цей елемент переходить із пружного в пластичний стан. Якщо матеріал зміцнюється при пластичному деформуванні, то відповідна межі пружності деформація ε е збільшується на величину D e е (Мал. 7):
(36)
Обчислення межі пружності по деформаціям , Досягнутого на кроці k визначається підсумовуванням:
. (37)
Мається на увазі, що в пружній області межа пружності не змінюється, його збільшення не обчислюються і дорівнюють нулю.
Накопичена пластична деформація визначається різницею інтенсивностей повної деформації e i та деформації e e, відповідної межі пружності:
(38)
Викладаються в подальшому ітераційні методи для досягнення задовільної збіжності вимагають дотримання безперервності і гладкості кривої зміцнення. Тому наприкінці пружного ділянки кривої зміцнення введений нелінійно пружний ділянку, на якій модуль зміцнення обчислюється за формулою:
, (39)
де - Інтенсивність деформації, відповідна межі пропорційності.
Співвідношення (39) висловлює пропорційну зміну модуля зміцнення при переході від пружного стану до пластичного. Межа пружності по напруженням в цьому випадку буде визначатися співвідношенням
, (40)
де e еp - деформація в області нелінійної пружності:
.
Вектор збільшень компонент тензора напруги на кроці k в пластичному стані визначається за збільшення компонент деформації:
. (41)
Вектор компонент напруги на кроці k у пружному та пластичному стані підсумовується за приростів:
(42).
Інтенсивність напружень визначається за компонентами тензора напруги s ij:
. (43)
Рис. 2. Зміна межі пружності за деформаціями при зміцненні
Якщо інтенсивність деформації зменшилася:
, (44)
то матеріал розвантажується і переходить в пружний стан. При порушенні нерівності (2.44) знову відбувається перехід елемента в пластичний стан.
5. Оцінка пошкоджуваності заготовок
Для оцінки деформованості і прогнозування руйнування заготовок в процесах обробки тиском отримала розвиток феноменологічна теорія руйнування, використання якої грунтується на отриманих досвідченим шляхом діаграмах пластичності та інформації про напружено-деформований стан в процесах обробки металів тиском.
Оцінку деформованості заготовок, а також розрахунок граничних технологічних параметрів проводять за допомогою деформаційних критеріїв, в основу яких покладені обмеження, що накладаються на деформації. При цьому для процесів, що супроводжуються монотонним, але складним деформуванням, в якості міри ушкоджень приймають зазвичай деяку скалярну характеристику.
Якщо впливом історії деформування знехтувати, то можна використовувати критерій Смирнова-Аляева:
(45)
Або, нормируя на одиницю, отримаємо міру пошкоджень y:
(46)
де e p (h) - гранична деформація в момент появи перших тріщин, які виявляються візуально; h - показник напруженого стану:
(47)
s - середнє нормальне напруження; s i - інтенсивність напружень.
Для обліку впливу історії деформування та використання співвідношення (47) для простого навантаження, приймемо за міру пошкоджень y вираз (критерій Колмогорова):
, (48)
де - Ступінь деформації до розглянутого моменту;
- Гранична деформація, обумовлена по діаграмах пластичності відповідних матеріалів.
Додавання в кінцево-елементну модель критерію деформованості дозволило проводити контроль на руйнування заготовки під час деформування, а також прогнозувати стан готового виробу.
6. Взаємодія заготовки з інструментом
Заготівля являє собою сукупність вузлів, пов'язаних між собою кінцево-елементної сіткою. Далі "вузлами" будуть називатися вузли звичайно-елементної сітки заготовки. Поведінка вузлів описується співвідношеннями МСЕ. Інструмент є абсолютно жорстким тілом, обмеженим непроникними для вузлів заготівлі кордонами. Межі інструменту апроксимуються прямолінійними і радіусними ділянками, зістикованими один з одним в єдину поверхню. Надалі під ім'ям "кордон" буде розумітися кордон інструменту. Інструмент може бути нерухомим і рухомим. Для адекватного опису технологічних процесів штампувального виробництва математична модель дозволяє використовувати один рухомий інструмент і кілька нерухомих (наприклад, матриця і оправлення).
Нерухомий інструмент зафіксований в просторі. У разі потрапляння вузла заготовки на кордон такого інструменту на його (вузла) переміщення накладаються обмеження - граничні умови, такі, що вузол може рухатися тільки вздовж кордону або від неї. Математично ця модель реалізується наступним чином.
У точці дотику вузла кордону обчислюється кут нахилу кордону. У осесиметричної задачі горизонталь паралельна радіальної осі, а вертикаль - осі симетрії тіла. Якщо межа в цій точці паралельна осьовому напрямку, то вузлу забороняються радіальні переміщення. Якщо ж межа паралельна радіальному напрямку, то вузлу забороняються осьові переміщення. У випадку з похилій кордоном встановлюється зв'язок між радіальної і осьової ступенями свободи вузла:
, (49)
де r - радіальна рівень свободи; z - осьова; k - тангенс кута нахилу кордону.
У матричному представленні це означає, що в матрицю жорсткості додається додатковий рядок з одиницею на позиції, що відповідає номеру забороненої ступеня свободи. Для похилій кордону використовується співвідношення:
(50)
Тобто в осередки, відповідні ступенями свободи r і z даного вузла, вносяться стоять при них в рівнянні (50) коефіцієнти. Таким чином, задовольняється рівняння зв'язку (49).
Рис. 3. Обмеження нерухомих кордонів
Рухома кордон переміщається в заданому напрямку (вертикальному або горизонтальному, в залежності від виду процесу) на величину D H на кожному кроці:
(51)
де H - повний хід інструменту; i - число кроків вирішення.
Якщо після чергового переміщення кордону будь-який вузол (або декілька) може виявитися "в тілі" інструмента, то йому призначається примусове переміщення на величину D h:
або
(52)
Тут R, Z - координати рухомий кордону; r, z - координати вузла.
Таким чином, переміщення вузла в напрямку перпендикулярної кордоні ступеня свободи буде забезпечувати його рух разом з кордоном, включаючи момент дотику. У цьому випадку переміщення вузла перпендикулярно кордоні не забороняється, а підпорядковується равенствам:
або
(53)
Для системи рівнянь це означає внесення у відповідні даному вузлу осередки матриці жорсткості коефіцієнтів 0 і 1 в залежності від напрямку руху і в праву частину величини D h.
Рис. 4 Обмеження рухливих кордонів
Рухома похила межа моделюється так само, як рухлива вертикальна або горизонтальна. Виникаюча при цьому помилка, пов'язана із зсувом вузла вздовж кордону мінімізується послідовними наближеннями (див. Рис. 5).
Рис. 5. Процедура уточнення положення вузла при зміщенні його рухомим інструментом
7. Тертя
При пластичному формозміни на межі контакту матеріалу та інструменту виникає сила тертя. При цьому напрямок течії матеріалу залежить від величини сили тертя, тобто напрямок сили тертя на межі контакту заздалегідь не відомо і може змінюватися.
Модуль напруги тертя визначається законом Кулона (54):
, (54)
де s N - нормальна напруга на кордоні інструменту; m - коефіцієнт тертя ковзання; t K - дотичне до кордону напругу.
Рис. 6. Контакт кінцевого елемента заготовки з інструментом
При пластичному плині дотичне напруження в елементах t K, що контактують із зовнішніми тілами (матриця або пуансон), не повинно перевищувати межу текучості при зсуві t S:
(55)
Для виконання цієї умови запропоновано наступний алгоритм розрахунку. Для елемента, що лежить на кордоні, визначається дотичне напруження t K, за яким обчислюється вузлова сила тертя:
(56)
де S K - площа контакту елемента з кордоном інструменту. Якщо до вузла примикають два елементи, що лежать на кордоні, то вузлова сила тертя виходить підсумовуванням сил, обчислених за формулою (56). Для визначення напряму сили тертя реалізується наступний алгоритм. Спочатку виконується крок навантаження без урахування тертя (F TP = 0). Потім за результатами цього кроку в кожному вузлі на межі контакту визначається нормальна вузлова сила на межі контакту, приріст переміщення вздовж неї і, описаним вище способом, сила тертя.
У процесі деформування вузли, що ковзають по кордоні інструменту, можуть зупинятися і потім змінювати напрямок руху. У нерухомому стані вузлова сила тертя може виявитися менше розрахованої за формулами (54) - (56). Для цього введено послідовне уточнення сили F TP додаванням до неї величини нев'язки D F TP:
(57)
Тут n - номер ітерації при уточненні сили тертя.
Початкове значення нев'язки визначається силою тертя:
(58)
Наступні значення нев'язки можуть міняти знак:
(59)
Ця формула забезпечує зростання сили тертя в напрямку протилежному руху. Якщо сила тертя зростає настільки, що змінює напрямок руху вузла, то ця ж формула забезпечить зменшення сили тертя.
Збіжність ітераційного процесу забезпечується зменшенням величини нев'язки при зміні напрямку руху:
(60)
Крім того, зростання сили тертя F TP обмежений її граничним значенням , Обчисленими за формулами (54) - (56). Описана итерационная процедура наближає силу тертя до значень, зображеним на графіку (Мал. 7).
Рис. 7. Залежність сили тертя від напрямку руху вузла
Основні висновки
1. Наведена математична модель, зроблена на базі методу скінченних елементів, є найбільш універсальною і адекватною з точки зору оцінки протікають в ній процесів.
2. Модель рухомий похилій кордону дозволяє більш точно уявити процес, а також оцінити картину деформування і течії матеріалу на радіусах заокруглення інструментів.
3. Створення на початковому етапі вирішення задачі трьох однакових копій матриці жорсткості і наступне в ході виконання накладення на дві з них умов, пов'язаних з переміщеннями вузлів, і використання третього матриці без зміни для визначення вузлових сил за переміщенням дозволило скоротити час обчислень, тому що створення копії є більш швидким процесом, ніж формування матриці жорсткості.
4. Введення в модель обліку пошкоджуваності заготовки дозволило не тільки проводити контроль за руйнацією заготовки безпосередньо під час протікання процесу, а й прогнозувати стан виробу, отриманого змодельованим методом.