Тут наводиться виклад спрощеного алгоритму вирішення плоскої задачі механіки деформованого твердого тіла методом кінцевих елементів у пакеті Mathcad, опублікованого в статті в журналі Exponenta.Pro (№3, 2003 р.), а також на форумі Exponenta.ru. Дату посту ставлю відповідну.

Розглядається найпростіший і водночас найпоширеніший варіант із розбиттям області на трикутні елементи. (Під час справи орієнтувався на алгоритм, викладений у книзі Фадєєв А.Б. Метод кінцевих елементів у геомеханіці. - М: Надра, 1987).

1. Підготовка вихідних даних.

Оскільки необхідно задати інформацію про кожен з елементів та вузлів розрахункової області, найзручніше використовувати для підготовки даних табличний редактор Excel, тим більше, що в Mathcad передбачена можливість імпорту даних із файлів формату.prn. Створюється в Excel два файли з таблицями, що містять відомості про елементи та вузли. Структура таблиць та розмірності величин у них наведено на рис. 1 і 2. У таблиці даних про вузли є два стовпці спеціальних змінних Рх і Ру, яким присвоюється ознака фіксації переміщення вздовж осі 0х або 0у відповідно (приймає значення 1, якщо задано нульове переміщення та 0 – якщо переміщення невідоме).

Мал. 1. Структура таблиці вихідних даних із інформацією про елементи.

Мал. 2. Структура таблиці вихідних даних з інформацією про вузли та задані вузлові сили та переміщення.

Для збереження таблиць у потрібному форматі, вибираємо Файл->Зберегти як…, вказуємо у відповідних полях ім'я файлу та тип – Форматований текст (розділювач – пробіл). Після натискання кнопки Зберегти, натиснути у вікні Так. Таким чином, отримуємо файли з іменами, наприклад EL_1.prn і KN_1.prn.

2. Завантаження даних у Mathcad. Підготовка змінних.

Для зручності нумерації елементів масивів далі в книзі Mathcad індекс перших елементів масивів встановлюється рівним одиниці:

Для отримання даних з файлів у Mathcad використовується функція READPRN(“filename.prn”) (можлива вказівка ​​повного шляху до файлу, інакше використовується поточна папка, шлях до якої можна дізнатися за допомогою функції CWD).

Припустимо, раніше створені файли знаходяться в папці DATA на диску D: Надається їх вміст матрицям DEL і DKN:

Надамо відповідним змінним значення з матриць:

Для перевірки правильності вихідних даних та використання в подальшому розрахунку необхідно сформувати вектор вузлових сил з урахуванням дії масових сил, вектор заданих переміщень, ознак фіксації переміщень і розрахувати площі елементів.

Площу n–го елемента зручно задати як функції користувача (у векторі V перераховуються глобальні номери вузлів елемента):

Загальна площа розрахункової галузі:

У масові сили перераховується вага елементів, порівну на кожен із їхніх вузлів:

Вузлові сили, переміщення та їх ознаки розміщуються у векторах послідовними парами значень: на парних позиціях вертикальні компоненти, на непарних горизонтальні:

3. Розрахунок матриці жорсткості системи.

Матриця жорсткості системи виходить шляхом об'єднання матриць жорсткості елементів [K] , які, в свою чергу, розраховуються за таким виразом.

Де Delta – площа елемента; [B] - матриця похідних функцій форми (функція впливу вузлів), [D] - матриця зв'язку напруг та деформацій:

Площа елемента обчислюється раніше заданою функцією користувача A(n). Матрицю [D] також зручно поставити як функції користувача; для умов плоскої деформації вона матиме вигляд

Матриця [B] пов'язує між собою переміщення вузлів елемента з його деформацією:

(вирази для функцій форми Nj, Nk виходять шляхом кругової підстановки індексів у порядку i, j, k),

i, j, k – номери вузлів елемента, xi, j, k, yi, j, k – координати вузлів.

Після нескладних перетворень матрицю [B] можна подати у вигляді

Матрицю [B] представимо як функції користувача, задавши попередньо допоміжну матрицю P , визначальну порядок перестановки індексів у функціях форми:

Матриця жорсткості системи обчислюється в наступному програмному блоці:

(Об'єднання матриць жорсткості елементів у МЖС проводиться за таким правилом: член МЖС Kci,j є сумою членів Ki,j з матриць жорсткості всіх елементів, що примикають до вузла з і-й ступенем свободи).

4. Розв'язання системи рівнянь

Після цього і-й стовпець та і-й рядок МЖС, а також і-й невідомий член у векторі сил можуть бути видалені. Для видалення рядків та стовпців із МЖС використовуємо субматриці, задані функціями користувача; М11 – видаляє перший рядок і стовпець, Mnn – останні, МI-IV – проміжні.

Таким чином, потрібно вирішити систему лінійних рівнянь алгебри (СЛАУ) . У разі можливості системи Mathcad дозволяють сильно спростити завдання. Для цього передбачено функцію lsolve(M,V) для знаходження вектора рішення СЛАУ, коефіцієнти якої містяться в масиві М, а вільні члени – у векторі V.

Програмний модуль ліворуч повертає «на місця» у загальному векторі задані вузлові переміщення, які раніше з нього видалені. Другий блок створює два вектори з осьовими компонентами вузлових переміщень.

5. Знаходження осьових деформацій та напруг в елементах

Знаючи отримані вузлові переміщення, можна розрахувати для кожного елемента деформації та напруги за співвідношеннями, згаданими в п.3 (сигма та епсілон):

У кожному елементі також підраховуються головні напруги та кут між віссю 0у та вектором максимальної головної напруги. Щоб уникнути поділу на нуль у рядку обчислення кута використано умовний вираз, який у разі рівності нулю виразу у знаменнику дробу надає куту значення .
.

6. Збереження результатів.

Розрахунок за вищевикладеною процедурою займає досить нетривалий час (наприклад, на ПК з процесором Pentium-IV-1300 MHz; 128 MB RAM час розрахунку для області зі 119 елементів 95 вузлів становить ~3 сек.), тим не менш, бажано зберегти результати для наступного аналізу.

Для цього сформуємо матриці, що характеризують ПДВ та поле переміщень, записавши в них також координати центрів елементів та вузлів:

(Для знаходження центрів елементів використана функція mean(), що повертає середнє значення елементів вектора)

Для запису даних у файл Mathcad передбачена функція WRITEPRN(«filename.prn»); перед її використанням можна задати попередньо кількість знаків після коми змінної PRNPRECISION та ширину стовпця у файлі змінної PRNCOLWIDTH:

Мал. 3. Розрахункова схема та її кінцево-елементне уявлення.

В даному випадку при розбитті на трикутні елементи вийшла мережа з 95 вузлів та 119 елементів. Нумерація – довільна.

Усі види навантаження, що діють на досліджувану область та формують у ній певний напружено-деформований стан, призводять до статично еквівалентних сил, прикладених у вузлах.

У силу симетрії граничні умови по переміщенням такі: горизонтальні компоненти уздовж вертикальної (x=0) і вертикальні уздовж горизонтальної (y=0) сторін квадрата дорівнюють нулю. Невідомі переміщення всіх вузлових точок всередині масиву, на контурі виробітку і на межі області.

Результати розрахунку можна представляти у вигляді епюр (рис. 4), ізоліній напруги або переміщень (рис. 5, а), поверхонь рівня (рис. 5, б). Збереження та подання результатів розрахунку у вигляді векторів (матриць) дозволяє зробити це без труднощів.

Мал. 4. Епюри напруги вздовж горизонтальної осі (для згладжування значень значення напруги призводять до центрів прямокутників, складених із сусідніх трикутників).

Мал. 5. Приклади візуалізації результатів розрахунку.

Література:

1. Фадєєв А.Б. Метод кінцевих елементів у геомеханіці. - М.: Надра, 1987. - 221с.
2. Єржанов Ж.С., Карімбаєв Т.Д. Метод кінцевих елементів у задачах механіки гірських порід. - Алма-Ата: Наука, 1975. - 239 с.
3. Зінкевич О. Метод кінцевих елементів у техніці: Пер. з англ. - М.: Світ, 1975. - 542 с.
4. Норрі Д., де Фріз Ж.. Введення в метод кінцевих елементів: Пер. з англ. - М.: Світ, 1981. - 304с.
5. Carlos A. Felippa. Introduction to Finite Element Methods. – Департамент Аерокосмічної архітектури і центру Аерокосмічними структурами University of Colorado, Boulder. - 2001.
6. Kyran D. Mish, Leonard R. Herrmann, LaDawn Haws. Finite Element Procedures in Applied Mechanics (траплялося десь в інеті).
7. Зенкевич О., Морган К. Кінцеві елементи та апроксимація. - М.: Світ, 1986. - 318 с.
8. Зенкевич О., Чанг І. Метод кінцевих елементів теорії споруд та в механіці суцільних середовищ. - М.: Надра, 1974. - 240 с.

Посилання:

• http://www.fea.ru/ ...Сайт FEA.RU присвячений актуальним проблемам кінцево-елементної механіки та комп'ютерного інжинірингу (CAE), МКЕ та розрахункам на міцність;
• http://www.cae.ru/ форум CAD та CAE-систем у тому числі теоретичні та прикладні аспекти КЕ моделювання та вирішення завдань механіки деформованого твердого тіла. Механіка конструкцій, машин, споруд та установок;
• - потужний каталог ресурсів, що стосуються МКЕ;
• http://www.isib.cnr.it/~secchi/EdMultifield/ - Сайт програми для кінцевоелементних розрахунків з непоганим описом методу.

Існує безліч програм, які вважають за МКЕ. Не вдаючись до подробиць того, чому метод такий гарний і широко застосовний, заглянемо в процес розрахунку зсередини. Здавалося б усе просто, чому не спробувати зібрати свій велосипед, тобто. зробити свою програму. На першому етапі налагоджувати, тестувати та налаштовувати розрахунок можна у MathCAD. Пізніше вже налагоджений алгоритм розрахунку для зручності введення даних та аналізу результатів можна буде переписати на C#, прикрутивши трохи графіки.

З чого починати? Так як моє глобальне завдання - моделювання ґрунту, то і починати розрахунки буду з задач теорії пружності.


Ось приклад завдання, яке необхідно розібрати. Пружні трикутні найпростіші КЕ. Схема складена та вирішена у програмі FEMmodels 2.0. Повторимо це в MathCAD.

1. Дискретизація області,
тобто розбиття цієї області на частини, визначення "вузлових" точок. З системи з нескінченним числом ступенів свободи складаємо систему з кінцевим числом вузлів і відповідно до ступенів свободи.

2. Визначення апроксимуючих функцій елемента.
Між вузлами значення функцій (у нашому випадку переміщення X і Y) змінюються за заданими нами законами, що апроксимують функціям.

3. Упорядкування рівнянь, що описують всю систему.
Як невідомі значення функцій у вузлах (у разі вийде система лінійних рівнянь СЛАУ).

4. Розв'язання рівнянь
та визначення вузлових значень та інших невідомих.

Я почну трохи не з розбиття області, а з другого пункту – визначення функцій кінцевого елемента. Найпростіший кінцевий елемент для розрахунку плоскої задачі теорії пружності - трикутник з лінійною функцією апроксимації:

Мал. 1. Апроксимуюча функція та отримання коефіцієнтів для неї.
Так як у нас два ступені свободи в кожному вузлі (X і Y), додається ще одна аналогічна функція.
Суть усіх маніпуляцій - отримання залежності між переміщеннями вузлів елемента і деформацій, що виникають у ньому. Так як у нас 6 компонентів переміщень і 3 деформацій, зв'язок здійснюється через деяку матрицю Bрозмірністю 3х6 (матриця похідних від функцій форм). Це перша матриця для побудови елемента.
Також потрібна матриця, що виражає взаємозв'язок деформацій та напруг (матриця D). Для пружного тіла ця залежність є узагальненим законом Гука.

Ще один маленький відступ від теми, матриця D для плоского напруженого стану іншого виду. Коли потрібно розрахувати, наприклад, основу насипу під залізничною колією або основу протяжної будівлі, тоді можна вважати плоску задачу, оскільки деформації вздовж насипу або будівлі можна прийняти нульовими. Для отримання D дорівнює e z =0. Якщо ж рахувати стіну будівлі, на яку діють сили тільки в площині стіни, теж можна вважати плоску задачу, тільки там деформації з площини перерізу будуть, а напруги немає, вважаємо sigma z =0.

Загальна матриця елемента K e:= B T D B V

Не переказуватиму математичні основи цього висновку, розповім короткий фізичний зміст.
Приклад матриці K e:

Число рядків і стовпців відповідає числу ступенів свободи. K i,j = зусилля у напрямі ступеня свободи j від додатку одиничного переміщення за напрямом ступеня свободи i. Тоді, наприклад, для нашого елемента в якості перевірки можна скласти парні/непарні елементи по якомусь рядку або стовпцю, за змістом нашого елемента це будуть реакції у закріпленнях Х або Y відповідно і їх сума природно дорівнює нулю. Матриця є виродженою, за фіз. сенс це означає, що незакріплений елемент має невизначені зусилля/реакції у вузлах.

Далі простіше. Потрібно зібрати глобальну матрицю системи із окремих елементів. За всіма ступенями свободи системи (рядки та стовпці матриці К) записуємо відповідні реакції від окремих елементів. Найдовше порався з цим перетворенням, а в результаті ось такий простий алгоритм з 5 вкладеними циклами:

Далі ще простіше, збираємо вектори сил, для всіх незакріплених ступенів свободи P, з Довикреслюємо рядки та стовпці із закріпленими ступенями свободи та отримуємо систему лінійних рівнянь: K * u = P; Вирішуємо u = K -1 P навіть не особливо замислюючись поки що про неекономічність цього методу в плані обчислювальної потужності, бо завдання мала.

Найнеприємнішим моментом рішення в маткаді залишилася незручність введення вихідних даних та аналіз результатів. Однак усі процедури алгоритмізуються і, наприклад, функція розміщення всіх закріплень займає 8 рядків, а 11 рядків складається список з n x n x 2 елементів (242шт у прикладі).

Мої два наступні завдання: елементи з більш складною апроксимацією, що дозволяють зменшити кількість елементів та уточнити рішення, і основне, нелінійні елементи. У цьому випадку матриця K буде залежною від переміщень і рішення суттєво ускладнюється. K(u)*u=P(u). У випадку вектор зовнішніх сил теж залежить від переміщень u.

Джерела знань:
1. Лекції 2008 р. на кафедрі ОІФ ПГУПС. Шашкін К.Г.
2. Сегерлінд "Застосування методу кінцевих елементів" (1979)
3. А.Л. Розін

Міністерство освіти та науки Російської Федерації

Державний освітній заклад

вищої професійної освіти

«Тихоокеанський державний університет»

Методичні вказівки та контрольні завдання до виконання

лабораторної роботи з курсу «Аналітичні та чисельні методи вирішення рівнянь математичної фізики» для студентів, які навчаються у магістратурі Хабаровськ Видавництво ТОГУ 2011 УДК 539.3/6. (076.5) Розв'язання двовимірної задачі теплопровідності методом кінцевих елементів у MATHCAD: методичні вказівки та контрольні завдання до виконання лабораторної роботи за курсом «Аналітичні та чисельні методи розв'язання рівнянь математичної фізики» для студентів, які навчаються в магістратурі/уклад. Л. М. Іванніков. - Хабаровськ: Вид-во Тихоокеан. держ. ун-ту, 2011.

Методичні вказівки складено на кафедрі «Механіка твердого тіла, що деформується». Включають зміст лабораторної роботи та рекомендації до вивчення розділів курсу «Аналітичні та чисельні методи вирішення рівнянь математичної фізики», необхідних для її виконання, список рекомендованої літератури та завдання для лабораторної роботи.

Друкується відповідно до рішень кафедри «Механіка деформованого твердого тіла» та методичної ради інституту будівництва та архітектури.

© Тихоокеанський державний університет,

ЗАГАЛЬНІ ПОЛОЖЕННЯ

Метою лабораторної є засвоєння алгоритму розрахунку двовимірних завдань теплопровідності методом кінцевих елементів.

РІВНЯННЯ ПЕРЕНОСУ ТЕПЛА

Рівняння плоскої задачі теплопровідності має вигляд 2T (x, y) 2T (x, y) Ky Q(x, y) 0, (1) Kx x 2 y кВт де K x, K y – коефіцієнти теплопровідності у напрямку осей x, y , ;

(м К) T (x, y) – потрібна функція температури; Q(x, y) – джерело тепла усередині тіла, кВт. Q(x, y) 0, якщо тепло підводиться до тіла.

м Граничні умови ставляться двох типів:

T TГ (Г), 1. (2) якщо температура T відома на деякій частині кордону Г, де TГ (Г) – відома температура в точках кордону, яка залежить від координат точок поверхні s на межі Г;

T (x, y) T (x, y) l Ky m h (T T) q (x, y) 0, 2. (3) K x x y якщо на частині поверхні Г 1 відбувається конвективний теплообмін, що характеризується величиною h (T T), або заданий потік тепла q(x, y) на частині поверхні Г 2, причому Г Г1 Г 2. Позначення в (2) та (3): h – коефіцієнт кВт теплообміну; T (x, y) - невідома температура на кордоні, К; T - (м2 К) відома температура навколишнього середовища, К; l, m - напрямні косинукВт си; q(x, y) – відомий потік тепла, вважається позитивним, якщо м тепло втрачається тілом. Потік тепла та конвективна тепловіддача на тому самому ділянці не можуть діяти одночасно.

Якщо має місце теплоізольована межа, то потік тепла дорівнює нулю і конвективний теплообмін відсутній, тоді гранична умова запишеться так:

dT 0, dn де n – зовнішня нормаль до межі області, що розглядається.

ФУНКЦІОНАЛ РІШЕННЯ ЗАВДАННЯ ТЕПЛОПРОВІДНОСТІ

Рішення рівняння (1) по області s з граничними умовами (2) і (3) на Г еквівалентно віднайденню мінімуму функціоналу При розв'язанні задачі МКЕ область s розбивається на n підобластей (кінцевих елементів), які зазвичай приймаються у формі трикутників (рис. 1) . Далі всі формули наводяться для трикутних КЕ. Функціонал записується як сума вкладів всіх кінцевих елементів області. Тоді (4) набуде вигляду провідності.

Або Представимо температуру, що змінюється в межах КЕ, через вузлові значення:

де [N(e)] - матриця функцій форми КЕ, що враховує розподіл температури в межах КЕ.

Тоді де [B(e)] - матриця градієнтів функцій форми КЕ.

Для кожного КЕ тепер можна записати вклад кожного КЕ у вираз для функціоналу (4):

Мінімум функціоналу (4) вимагає виконання наступної умови:

Для окремого КЕ отримаємо де матриця теплопровідності КЕ має вигляд а вектор зовнішнього впливу буде Для всієї області отримаємо або де Рівняння (6) є основним рівнянням для вирішення задачі теплопровідності методом кінцевих елементів.

ДВОМІРНИЙ СИМПЛЕКСЕЛЕМЕНТ

Для вирішення плоского завдання теплопровідності використовується трикутний КЕ з прямолінійними сторонами (див. рис. 1). Нумерація вузлів проводиться проти годинникової стрілки, починаючи з деякого вузла, що позначається одиницею.

Нумерація сторін КЕ наведено на рис. 1.

Вузлові значення температури позначаються T1, T2, T3. Температура у точці КЕ з координатами x, y визначається за формулою Нижче наводяться функції форми, що застосовуються для цього КЕ.

Площа КЕ обчислюється за відомою формулою Коефіцієнти, що входять до функції форми, залежать від координат вузлів, наведені нижче:

ЗАСТОСУВАННЯ ЧОТИРИКУТНИХ КЕ ДЛЯ ГЕНЕРАЦІЇ сітки

Для попереднього нанесення сітки з великим осередком (розбивкою області на зони) використовуються чотирикутні квадратичні елементи (рис. 2).

На кожній стороні КЕ вводиться три вузли.

На рис. 2 показані місцеві відносні координатні осі, у яких вузол 7 (1, 1). Нумерація вузлів такого КЕ починаючи з вузла 1 проводиться проти годинникової стрілки. Вузли 2, 4, 6, 8 можуть розташовуватися в довільній точці відповідної сторони, що дозволяє будувати більш густу сітку поблизу точкових впливів. Надалі кожна сторона такого КЕ розбивається на задану кількість ділянок. Нумерація вузлів проводиться так: по вертикалі від вузла з координатами (1, 1) вниз по осі і зліва направо по осі. Таким чином, великі елементи діляться на дрібніші, які, у свою чергу, меншою за довжиною діагоналлю розбиваються на трикутні КЕ. Трикутні ділянки зони також представляються як чотирикутних квадратичних елементів (рис. 3).

Мал. 3. Подання трикутної області у вигляді чотирикутного квадратичного елемента

МАТРИЦЯ ТЕПЛОПРОВІДНОСТІ КЕ

Для трикутного КЕ матриця теплопровідності має вигляд де L12, L23, L31 – довжини відповідних сторін КЕ. Останні три члени враховують конвективний теплообмін з кожної сторони КЕ. Так як КЕ входить складовою в розглянуту область, то конвективний теплообмін зазвичай відбувається по одній або двох сторонах КЕ.

ВЕКТОР ЗОВНІШНІХ ВПЛИВ НА КЕ

Зовнішніми (відомими) впливами є:

1. Джерело тепла всередині КЕ постійної інтенсивності Q(e).

2. Приплив тепла за рахунок теплового потоку q(e).

3. Конвективний теплообмін не більше ніж по обидва боки КЕ з коефіцієнтом теплообміну h(e).

4. Точкове джерело тепла Q * (X 0, Y0), що знаходиться всередині КЕ.

Вектор зовнішніх впливів на КЕ має вигляд

ГРАДІЄНТИ ТЕМПЕРАТУР І СЕРЕДНЯ ТЕМПЕРАТУРА ПО КЕ

Градієнти температур та середня температура по КЕ обчислюються за такими формулами:

ПОРЯДОК РІШЕННЯ ЗАВДАННЯ ТЕПЛОПРОВІДНОСТІ

НАНЕСЕННЯ сітки вузлів на розглянуту область

Область розв'язання задачі міститься в системі глобальних координат X, Y. Розглянута область повинна бути покрита сіткою вузлів. Чим менше буде комірка сітки, тим точніше буде вирішення задачі. Нанесення сітки проводиться згідно з 2 етапами.

І етап. Розглянута область розбивається на ряд прямокутних та трикутних зон (чотирикутні квадратичні елементи). Зони нумеруються у довільному порядку. Для кожної такої зони задаються 8 вузлових точок (по три на кожній стороні, включаючи кутові точки). Для трикутної зони одна із сторін відповідає двом сторонам прямокутника (5 точок).

Таким чином, при розбивці на зони використовують чотирикутні квадратичні елементи.

Складаються такі таблиці вихідних даних:

а) Табл. 1 з'єднання зон, що визначає, які сторони зон контактують між собою.

З'єднання зон у розглянутій області. Таблиця 1.

У наведеній таблиці. 1 показано, що зона 1 контактує тільки з зоною по першій стороні, зона 2 контактує з зоною 1 по першій стороні та з зоною 3 по четвертій стороні. Зона 3 контактує лише із зоною 2 з другого боку (рис. 4). Нумерація сторін залежить від орієнтації місцевих осей у відносних координатах, які показані малюнку жирними цифрами. На рис. 4 показано напрямок нумерації вузлів зон від початкового вузла Н.

Мал. 4. Формування таблиці з'єднання зон б). Табл. 2 координат вузлів, нанесених на межі зон, у прийнятій глобальній системі координат.

Координати вузлів на межі зон Таблиця 2.

в). Табл. 3, в якій вказується число смуг по вертикалі та горизонталі, на які розбивається кожна зона для отримання сітки з осередками менших розмірів.

Формування сітки з меншими за розміром осередками Таблиця 3.

Зона 1 розбивається на п'ять смуг по висоті та шість смуг по ширині.

г). Табл. 4, у якій кожної зони вказуються раніше нанесені вузли.

Номери вузлів попередньої сітки для кожної зони Таблиця 4.

У табл. 4 зазначено, що вісім вузлом другої зони мають такі номери при обході зони, що розглядається, проти годинникової стрілки.

ІІ етап. Далі в Mathcad реалізована програма "grid", в якій задається число смуг за висотою та шириною для кожної зони, що дозволяє розбити кожну зону на прямокутники набагато менших розмірів. Потім кожен з цих малих прямокутників меншою по довжині діагоналлю ділиться на два трикутники і вся область, що розглядається, покривається сіткою з трикутною коміркою.

В результаті роботи цієї програми видаються такі дані:

a). Число трикутних КЕ (Kol_Elem).

б) Наступні таблиці. 5, 6, 7.

Нумерація вузлів сітки по сторонах зон Таблиця Таблиця видається у формі матриці розміром (кількість смуг зони по висоті число смуг зони по ширині) для кожної зони, що спрощує побудову сітки.

Наведена матриця показує, що в зоні 3 на стороні 1 розташовуються вузли 23, 24, 25, 26; на боці 2 розташовуються вузли 26, 22, 1; на боці 3 – вузли 1, 16, 13, 10; на стороні 4 вузли 10, 19, 23. Обхід зони проти годинникової стрілки. Ця нумерація показана у наведеному нижче прикладі.

Таблиця Також можуть бути виведені таблиці, що зв'язують номер зони, номер КЕ і координати вузлів КЕ.

На схему області, що розглядається, вручну наноситься сітка з нумерацією КЕ та їх вузлів.

ФОРМУВАННЯ ВЕКТОРА ЗОВНІШНІХ ВПЛИВ

На підставі побудованої сітки для цієї області відзначаються:

а) Номери сторін, якими відбувається конвективний обмін тепла.

б) Номери вузлів, у яких задана температура.

в) Номери КЕ, в яких на їх сторонах, вузлах або всередині розміщуються зосереджені теплові джерела.

Складаються такі табл. 8, 9, 10.

Сторони області з конвективною тепловіддачею Таблиця Передбачається, що конвективний теплообмін можливий лише з обох боків КЕ із трьох.

Таблиця точкових джерел тепла Таблиця Таблиця величин температури у вузлах КЕ.

Таблиця градієнтів температур Gradx, Grady по осях Х та Y відповідно.

Таблиця середньої температури Тsred по кожному КЕ.

Розподіл температур по області, що розглядається, із зазначенням величин ізотерм.

ПРИКЛАД ВИКОНАННЯ ЛАБОРАТОРНОЇ РОБОТИ

У теплопровідному середовищі проходять 4 кабелі, як показано на рис. 5. Середовище має коефіцієнти теплопровідності K x K y 10. Коефіцієнт тепсм К лообміну на поверхні середовища h 5. По бокових сторонах розглянемо 2 К ванна середовище обмежена товстим шаром ізоляції. Температура повітря на поверхні середовища T300C. Температура нижнього шару середовища T200C.

Потужність випромінювання тепла кожним кабелем становить Q 200 Вт.

Потрібно:

ВКАЗІВКИ:

а) при виконанні лабораторної роботи врахувати симетрію області та симетрію температурного впливу;

б) розбити частину області, що розраховується, на три або чотири зони;

в) кожну зону розбивати від трьох до п'яти смуг за висотою та шириною для спрощення нанесення сітки на область.

РІШЕННЯ ЗАВДАННЯ

Враховуючи симетрію аналізованої області, в розрахунку будемо враховувати лише половину цієї області (рис. 6).

Помістимо розглянуту область у систему глобальних осей X і Y і розіб'ємо її на три зони, на сторони яких нанесемо вузли, вважаючи зони чотирикутними квадратичними елементами рис. 7. Пронумеруємо зони та вузли, обходячи область проти годинникової стрілки. Для визначення номерів сторін зон для кожної зони встановлюється система місцевих осей.

Мал. 7. Попередня розбивка області на зони Для більш точного вирішення завдання необхідно вузли на межі зон розташовувати ближче до точкових джерел тепла.

Складаються вихідні дані для призначених зон та вузлів (табл. 1, 2, 3, 4). Програма розрахунку видає табл. 5, 6, 7, що надають повну інформацію про трикутну сітку, нанесену на область, що використовується в подальшому розрахунку. За цими таблицями на аркуші будується сітка (рис. 8).

Мал. 8. Трикутна сітка, нанесена на область За отриманою сіткою проводиться облік зовнішнього температурного впливу та складаються табл. 8, 9, 10. Після чого у табличній формі виводяться результати розв'язання задачі та їх графічне подання на рис. 9 та

1. Р О З У Л Ь Т А ТИ РІШЕННЯ З СТВОРЕННЯ сітки КЕ

ВУЗЛИ сітки по кордонах ЗОН

ТАБЛИЦЯ КЕ

Координати вузлів КЕ

ФОРМУВАННЯ ВЕКТОРА

ЗОВНІШНІХ ВПЛИВ

stor

2. Р О З У Л Ь Т А ТИ РІШЕННЯ ЗАВДАННЯ

Градієнти температур та середня температура по КЕ перетури по області

ВАРІАНТИ ЗАВДАНЬ ДЛЯ ЛАБОРАТОРНОЇ РОБОТИ

У теплопровідному середовищі, як показано на схемі, проходять кабелі, що випромінюють тепло. Середовище має коефіцієнти теплопровідності K x та K y. Коефіцієнт теплообміну поверхні середовища h. На деяких ділянках середовище, що розглядається, обмежена товстим шаром ізоляції. Температура повітря окремих ділянках середовища, де відбувається конвективний теплообмін, Т. На деяких ділянках середовища задана температура Т.

Потужність випромінювання тепла кожним кабелем становить Q.

Потрібно, використовуючи вихідні дані для свого варіанту та схему завдання (табл. 11, рис. 11):

1. Визначити розподіл температури у заданій області.

2. Визначити градієнти температур та середню температуру по області.

3. Побудувати графіки зміни отриманих величин.

ВИХІДНІ ДАНІ

Вихідні дані до лабораторної роботи за варіантами Таблиця Ноh, анта Рис. 11. Схеми варіантів задній для лабораторної роботи

КОНТРОЛЬНІ ПИТАННЯ

Запишіть рівняння теплопровідності для двовимірного завдання.

Запишіть граничні умови для двовимірного завдання теплопровідності.

Запишіть повний функціонал розв'язання задачі теплопровідності.

Виведіть основне рівняння для вирішення двовимірної задачі теплопровідності шляхом кінцевих елементів.

5. Які кінцеві елементи використовуються для вирішення двовимірної задачі теплопровідності?

6. Як визначаються функції форми для двовимірного симплекселемента?

7. З якою метою використовуються чотирикутні квадратичні елементи?

8. Як вибирається система місцевих координат та проводиться нумерація сторін чотирикутного квадратичного елемента?

9. Запишіть матрицю теплопровідності для трикутного КЕ.

10. Як формується матриця теплопровідності для аналізованої області?

11. Як формується вектор зовнішніх теплових впливів для КЕ?

12. Як формується вектор зовнішніх впливів для цієї області?

13. Як визначаються градієнти температур та середня температура по КЕ?

14. Як проводиться нанесення сітки на область, що розглядається?

15. Які вихідні дані потрібно підготувати для нанесення сітки?

16. Які вихідні дані використовуються для побудови сітки та як вона наноситься на область?

17. Які дані необхідно внести на формування вектора зовнішніх теплових впливів?

18. Як зважити на знак величини точкового джерела тепла? Приплив тепла?

19. Які вихідні дані виходять у результаті розв'язання задачі теплопровідності?

1. Зенкевич О. Метод кінцевих елементів у техніці / О. Зенкевич. - М.:

Мир, 1975. - 452 с.

2. Сегерлінд Л. Застосування методу кінцевих елементів/Л. Сегерлінд. - М.:

Мир, 1979. - 392 с.

ЗАГАЛЬНІ ПОЛОЖЕННЯ ………………………………………………. Рівняння перенесення тепла……………….…………………………… … Функціонал розв'язання задачі теплопровідності……………………..... Двовимірний симплекс елемент…………………. .………………..…. Застосування чотирикутних КЕ для генерації сітки....……….… Матриця теплопровідності КЕ…………………….………………... Вектор зовнішніх впливів на КЕ….………………… ………… …... Градієнти температур і середня температура по КЭ…………… …… Порядок вирішення задачі теплопровідності в Mathcad 14…..….… ПРИКЛАД ВИКОНАННЯ ЛАБОРАТОРНОЇ РОБОТИ…….… ….. РІШЕННЯ ЗАДАЧІ…… ………………………………………………. Роздрукування розв'язання задачі……………………………………… …..... ВАРІАНТИ ЗАВДАНЬ ДЛЯ ЛАБОРАТОРНОЇ РОБОТИ…..... Контрольні питання………………………………… ………… ……. Бібліографічний список……………………………………… ……

РІШЕННЯ ДВОМІРНОЇ ЗАВДАННЯ ТЕПЛОПРОВІДНОСТІ

МЕТОДОМ КІНЦЕВИХ ЕЛЕМЕНТІВ У MATHCAD

Методичні вказівки та контрольні завдання для виконання лабораторної роботи з курсу «Аналітичні та чисельні методи вирішення рівнянь математичної фізики» для студентів, які навчаються у магістратурі.

Головний редактор А. А. Суевалова Редактор Т. Ф. Шейкіна Оператор комп'ютерної верстки Л. М. Іванніков Підписано до друку Папір писчий. Гарнітура "Таймс". Друк цифровий.

Ум. піч. л. Тираж 50 екз. Видавництво Тихоокеанського державного університету.

680035, Хабаровськ, вул. Тихоокеанська, 136.

Відділ оперативної поліграфії видавництва Тихоокеанського державного університету. 680035, Хабаровськ, вул. Тихоокеанська, 136.

Схожі роботи:

В.Б. Пономарьов А.Є. Замураєв АСПИРАЦІЯ ТА ОЧИЩЕННЯ ПРОМИСЛОВИХ ВИКИДІВ І СКИДАНЬ Федеральне агентство з освіти ГОУ ВПО Уральський державний технічний університет-УПІ В.Б. Пономарьов А.Є. Замураєв АСПИРАЦІЯ ТА ОЧИЩЕННЯ ПРОМИСЛОВИХ ВИКИДІВ І КИДІВ МЕТОДИЧНІ ВКАЗІВКИ ПО КУРСУ МАШИНИ ТА АГРЕГАТИ ПІДПРИЄМСТВ БУДІВЕЛЬНИХ МАТЕРІАЛІВ Науковий редактор – проф. техн. наук В.Я.Дзюзер Єкатеринбург УДК 666.9.001.575 (042.4) ББК 35.41в П Рецензенти: Пономарьов В.Б. П56 Аспірація та...»

«Міністерство освіти і науки Російської Федерації Федеральна державна бюджетна освітня установа вищої професійної освіти Тихоокеанський державний університет ВИРОБНИЦТВО КАМ'ЯНИХ РОБІТ програм підготовки спеціаліста...»

«МІНІСТЕРСТВО ОСВІТИ І НАУКИ РОСІЙСЬКОЇ ФЕДЕРАЦІЇ Федеральне агентство з освіти Державна освітня установа вищої професійної освіти РОСТІВСЬКИЙ ДЕРЖАВНИЙ БУДІВЕЛЬНИЙ УНІВЕРСИТЕТ ЗАТВЕРДЖЕНО на засіданні кафедри економіки та управління в будівництві2. МЕТОДИЧНІ ВКАЗІВКИ щодо виконання науково-дослідної роботи для студентів, магістрантів та аспірантів економічних спеціальностей та напрямків Ростов-на-Дону, УДК ​​69.003(07)...»

«База нормативної документації: www.complexdoc.ru Довідкові матеріали Матеріали та технології століття Добромислов А.Я., Санкова Н.В. Пластмасові труби та сучасні технології для будівництва та ремонту трубопроводів ПРОЕКТУВАННЯ, МОНТАЖ ТА ЕКСПЛУАТАЦІЯ СИСТЕМ КАНАЛІЗАЦІЇ З ПЛАСТМАСОВИХ ТРУБ ДЛЯ БУДІВЕЛЬ І МІКРОРАЙОНІВ РЕКОМЕНДАЦІЇ Глава 1004. З ПОЛІМЕРНИХ МАТЕРІАЛІВ ДЛЯ МОНТАЖУ ТРУБОПРОВІДІВ СИСТЕМ КАНАЛІЗАЦІЇ БУДІВЕЛЬ І МІКРОРАЙОНІВ 1.1. Внутрішня...»

«МІНІСТЕРСТВО ОСВІТИ РОСІЙСЬКОЇ ФЕДЕРАЦІЇ Державна освітня установа вищої професійної освіти Оренбурзький державний університет Кафедра технології будівельних матеріалів та виробів Т.І. ШЕВЦОВА МАТЕРІАЛОВЕДЕННЯ МЕТОДИЧНІ ВКАЗІВКИ З ВИВЧЕННЯ ДИСЦИПЛІНИ “МАТЕРІАЛОВЕДЕННЯ” Рекомендовано до видання Редакційно-видавничою радою Державної освітньої установи вищої професійної освіти “Оренбурзький державний університет...

«Міністерство освіти і науки Російської Федерації Санкт-Петербурзький державний архітектурно-будівельний університет Г. П. КОМІНА, А. О. ПРОШУТИНСЬКИЙ ГІДРАВЛІЧНИЙ РОЗРАХУНОК І ПРОЕКТУВАННЯ ГАЗОПРОВІДІВ Навчальний посібник Санкт-Петербург 2010 УДК 682.6. техн. наук, доц. М. А. Кочергін, головний спеціаліст відділу технічного нагляду Управління капітального будівництва ВАТ Газпромрегіонгаз; А. Г. Матвєєв, заст. генерального директора Інституту...»

«ДЕПАРТАМЕНТ ПРАЦІ ТА СОЦІАЛЬНОГО ПІДТРИМКИ НАСЕЛЕННЯ ЯРОСЛАВСЬКОЇ ОБЛАСТІ Реалізація обласної цільової програми Доступне середовище. Організація роботи органів соціального захисту населення та установ соціального обслуговування населення Ярославської області щодо соціальної реабілітації інвалідів ЗБІРКА ІНФОРМАЦІЙНИХ І МЕТОДИЧНИХ МАТЕРІАЛІВ Ярославль 2011 Реалізація обласної цільової програми Доступне середовище. Організація роботи органів соціального захисту населення та установ соціального обслуговування...»

«МІНІСТЕРСТВО ОСВІТИ І НАУКИ РОСІЙСЬКОЇ ФЕДЕРАЦІЇ Федеральне державне бюджетне навчальний заклад вищої професійної освіти АНІЯ З ПРАКТИЧНИХ ЗАНЯТТІВ І САМОСТІЙНОЇ ПРАЦІ СТУДЕНТІВ для направлення 081100.62 Державне та муніципальне управління очної та заочної форми навчання Тюмень,...»

«Федеральне агентство з освіти КАЗАНСЬКИЙ ДЕРЖАВНИЙ АРХІТЕКТУРНО-БУДІВЕЛЬНИЙ УНІВЕРСИТЕТ Н.С. ББК 24 До 78 Громаков М. З. Хімічний зв'язок: Методичні вказівки з хімії для студентів денної, заочної та дистанційної форм навчання, Казань: КДАСУ, 2007. -37с. Методичні вказівки містять основний інформаційний матеріал,...»

«МІНІСТЕРСТВО ОСВІТИ ТА НАУКИ УКРАЇНИ ХАРКІВСЬКА НАЦІОНАЛЬНА АКАДЕМІЯ МІСЬКОГО ГОСПОДАРСТВА В.І. ОСПІЩІВ ОСНОВИ МАРКЕТИНГУ Навчальний посібник (для студентів спеціальності 6.070101 – Транспортні технології) Харків Видавництво “ФОРТ” 2009 УДК 339.138(075.8) ББК 65.290-2я7 О75 Рецензенти: Іванілов, д.е.н., професор, зав. кафедрою економіки Харківського державного технічного університету будівництва та архітектури; Г.В. Ковалевський, д.е.н., професор кафедри маркетингу та ме...»

«Федеральна агенція з освіти Державна освітня установа вищої професійної освіти Володимирський державний університет Кафедра будівельних конструкцій МЕТОДИЧНІ ВКАЗІВКИ ДО ВИВЧЕННЯ РОЗДІЛУ ЗАЛІЗОБЕТОННІ КОНСТРУКЦІЇ КУРСУ КОНСТРУКЦІЇ МІСЬКИХ ЗДАНЬ. МИХАЙЛОВ В.І. ВОРОНІВ Володимир 2009 УДК 624.012.3/4 ББК 38.53 М54 Рецензент Кандидат технічних наук, доцент зав. кафедрою будівельних конструкцій Володимирського...»

«МІНІСТЕРСТВО ОСВІТИ РОСІЙСЬКОЇ ФЕДЕРАЦІЇ Державний навчальний заклад вищої професійної освіти – Оренбурзький державний університет О.С. КІЛОВ ОБРОБКА МАТЕРІАЛІВ ТИСКОМ У ПРОМИСЛОВОСТІ Рекомендовано Вченою радою державної освітньої установи вищої професійної освіти Оренбурзький державний університет як навчальний посібник для студентів, які навчаються за програмами вищої професійної освіти з...»

«Федеральна агенція з освіти Ульянівський державний технічний університет Хімія води Навчальний посібник для студентів нехімічних спеціальностей Укладачі: Л.В. Петрова, Є.М. Калюкова Ульяновськ 2004 УДК 541.1(075.8) ББК 24 Я7 Х 46 Рецензенти: Зав. науково-виробничим відділом ТОВ Трансбудкомплект, канд. техн. наук. І.А. Дорофєєв, Зав. кафедрою Хімія УлДПУ, доцент, канд. хім. наук І. Т. Гусєва Затверджено редакційно-видавничою радою університету як навчальний...»

«Міністерство освіти та науки Державна освітня установа вищої професійної освіти Пермський державний університет ПРИРОДНО НАУКОВИЙ ІНСТИТУТ Н.Г. Максимович, Є.А. Хайруліна ГЕОХІМІЧНІ БАР'ЄРИ ТА ОХОРОНА НАВКОЛИШНЬОГО СЕРЕДОВИЩА Навчальний посібник Перм 2011 УДК 504.06:550.4 ББК 20.18:26.30 M 18 Максимович, Н.Г. М18 Геохімічні бар'єри та охорона навколишнього середовища: навч. посібник/Н.Г. Максимович, Є.А. Хайруліна; Пермь. держ. ун-т. - Перм, 2011. - 248с.: іл. ISBN...»

«МІНІСТЕРСТВО ОСВІТИ ТА НАУКИ УКРАЇНИ ХАРКІВСЬКИЙ НАЦІОНАЛЬНИЙ УНІВЕРСИТЕТ МІСЬКОГО ГОСПОДАРСТВА імені О. М. БЕКЕТОВА ПОКРІВЕЛЬНІ ТА ГІДРОІЗОЛЯЦІЙНІ НАВЧАЛЬНЕ ПІДСУМКИ будівництво За редакцією В. Д. Жван Харків ХНУГГ 2013 УДК (075 ) ББК 38.654.3я73-6+38.637я73-6 К83 Автори: Жван Віктор Денисович – кандидат технічних наук, професор, професор кафедри Технологія...»

« МЕНЕДЖМЕНТ Навчальний посібник Для студентів вузів У двох частинах Частина 1 Кемерово 2008 2 УДК 65.018 (075) ББК 30.607я7 М 31 Рецензенти: Є.Г. Ягупа, канд. екон. наук, доцент, зав. кафедрою Економічна теорія та економіка підприємств КДСГІ; С.М. Бугрова, канд. екон. наук, доцент кафедри Економіка та організація машинобудівної...»

«Міністерство освіти та науки Російської Федерації ГОУ ВПО Тамбовський державний технічний університет О.В. УМНОВА, О.В. ЄВДОКИМЦІВ СТАЛЬНИЙ КАРКАС БУДИНКУ ПАВІЛЬЙОННОГО ТИПУ Затверджено Вченою радою університету як навчальний посібник Тамбов Видавництво ТДТУ 2008 УДК 624.014.2(075) ББК Н549 У545 Р е ц. Леденєв Генеральний директор ТОВ ФСК Тамбоврегіонбуд В.І. Скрильов Умнова, О.В. У545 Сталевий каркас будівлі павільйонної...»

«Міністерство освіти та науки Російської Федерації ФДБОУ ВПО Казанський державний архітектурно-будівельний університет Кафедра Економіки та підприємництва у будівництві Методичні вказівки Організація презентацій. Використання мультимедіа-ресурсів для студентів наступних спеціальностей та напрямів підготовки: 070603 Мистецтво інтер'єру, 190702 Організація та безпека руху, 190205 Підйомно-транспортні, будівельні, дорожні машини та обладнання, 270100.62...»

«МИНИСТЕРСТВО ОБРАЗОВАНИЯ И НАУКИ РОССИЙСКОЙ ФЕДЕРАЦИИ ФЕДЕРАЛЬНОЕ АГЕНТСТВО ПО ОБРАЗОВАНИЮ САНКТ-ПЕТЕРБУРГСКИЙ НАЦИОНАЛЬНЫЙ ИССЛЕДОВАТЕЛЬСКИЙ УНИВЕРСИТЕТ ИНФОРМАЦИОННЫХ ТЕХНОЛОГИЙ, МЕХАНИКИ И ОПТИКИ А.А. Воробйова НАВЧАЛЬНИЙ ПОСІБНИК ПО КУРСУ ГЕОІНФОРМАЦІЙНІ СИСТЕМИ ТЕРИТОРІАЛЬНОГО УПРАВЛІННЯ Санкт-Петербург 2012 1 Навчальний посібник присвячений геопросторовому моделюванню об'єктів за допомогою ГІС та використання семантів, що їх супроводжують. Крім того питанням збору та підготовки...»

«МІНІСТЕРСТВО ОСВІТИ ТА НАУКИ РОСІЙСЬКОЇ ФЕДЕРАЦІЇ КАЗАНСЬКИЙ ДЕРЖАВНИЙ АРХІТЕКТУРНО-БУДІВЕЛЬНИЙ УНІВЕРСИТЕТ Кафедра виробничої безпеки та права БЕЗПЕКА для студентів00 .65 та профілю 270804.62 Виробництво та застосування будівельних матеріалів, виробів та конструкцій Казань 2013 УДК 69.05 : 658.382 ББК До 66 До 66 Безпека життєдіяльності:...»

А.В. Ігнатьєв, Н.А. Михайлова, Т.В. Ерещенко МЕТОД КІНЦЕВИХ ЕЛЕМЕНТІВ І ЙОГО РЕАЛІЗАЦІЯ У СЕРЕДОВИЩІ MATHCAD Лабораторний практикум з дисципліни «Аналітичні та чисельні методи вирішення рівнянь математичної фізики» Волгоград 2010 Міністерство освіти і науки Російської Федерації Волгоградський державний технік. Ігнатьєв, Н.А. Михайлова, Т.В. Ерещенко МЕТОД КІНЦЕВИХ ЕЛЕМЕНТІВ ТА ЙОГО РЕАЛІЗАЦІЯ У СЕРЕДОВИЩІ MATHCAD Лабораторний практикум з дисципліни «Аналітичні та чисельні методи розв'язання рівнянь математичної фізики» Волгоград 2010 УДК 624.04:012.2.2. 3-018.2я73 І 266 Рецензія н ти: кандидат технічних наук М.М. Степанов, професор кафедри прикладної математики та обчислювальної техніки ВолгДАСУ; доктор технічних наук Н.Г. Бандурін, професор кафедри будівельної механіки ВолгДАСУ Затверджено редакційно-видавничою радою університету як навчально-практичний посібник Ігнатьєв О.В. І 266 Метод кінцевих елементів та його реалізація у середовищі Mathcad: лабораторний практикум з дисципліни «Аналітичні та чисельні методи вирішення рівнянь математичної фізики»/О.В. Ігнатьєв, Н.А. Михайлова, Т.В. Єрещенко; Волгогр. держ. архіт.будує. ун-т. Волгоград: ВолгДАСУ, 2010. 31 с. ISBN 978-5-98276-372-3 Міститься короткі теоретичні відомості, необхідні для виконання лабораторної роботи з дисципліни «Аналітичні та чисельні методи вирішення рівнянь математичної фізики», наведені варіанти індивідуальних завдань, приклади виконання завдань лабораторної роботи, а також досліджуваної теми. Для магістрів спеціальності АТ, СМ та ВІВ денної форми навчання. УДК 624.04:004.92(076.5) ББК 38.112я73+32.973-018.2я73 ISBN 978-5-98276-372-3 Державна освітня установа вищої професійної освіти «Волгоградський державний02, Волгоградський державний02 МЕТОД КІНЦЕВИХ ЕЛЕМЕНТІВ ТА ЙОГО РЕАЛІЗАЦІЯ У СИСТЕМІ MATHCAD Мета роботи: вивчити метод кінцевих елементів та отримати навички його реалізації в системі Mathcad. Основні поняття і концепція МКЕ Основна ідея методу Основна ідея методу полягає у поданні конструкції, що розраховується, у вигляді сукупності елементів простої форми, з'єднаних між собою в окремих точках. По суті, суцільне середовище з нескінченним числом ступенів свободи замінюється набором підобластей, що мають кінцеве число ступенів свободи. При такому підході шукані безперервні величини (переміщення, напруги, деформації тощо) всередині кожного кінцевого елемента (КЕ) виражаються за допомогою апроксимуючих функцій через вузлові значення цих величин. Розподілені зовнішні навантаження замінюються еквівалентними вузловими силами. У математичному плані завдання полягає у приведенні диференціальних рівнянь або енергетичного функціоналу, що описують конструкцію, що розглядається, до системи алгебраїчних рівнянь, вирішення якої дає значення шуканих вузлових невідомих. Метод кінцевих елементів має дуже широке поширення у практиці розрахунків на міцність, стійкість та коливання будівельних, машинобудівних, авіаційних конструкцій. За допомогою МКЕ можна успішно виконати аналіз широкого класу стрижневих систем (ферми, рами тощо), тонкостінних просторових конструкцій (плити перекриттів, оболонки покриттів тощо), масивних тривимірних тіл, а також комбінованих систем, що складаються з одновимірних. , двовимірних та тривимірних елементів. МКЕ відрізняє широка область застосування, інваріантність по відношенню до геометрії конструкції та фізичних характеристик матеріалів, відносна простота обліку взаємодії конструкцій з навколишнім середовищем (механічні, температурні, корозійні впливи, граничні умови тощо), високий ступінь пристосовності до автоматизації всіх етапів розрахунку . Метод має просту фізичну інтерпретацію і тісно пов'язаний із методом переміщень, який широко використовується у будівельній механіці. 3 На основі кінцево-елементного підходу розроблено велику кількість потужних програмних комплексів. Серед них можна відзначити такі як ABAQUS, ADINA, ANSYS, COSMOS/М, MSC/NASTRAN, ЛІРА, SCAD, STARK, СТАДІО. Більшість з них має велику бібліотеку кінцевих елементів і дає можливість виконувати розрахунки на міцність, стійкість і коливання, враховувати фізичну та геометричну нелінійності, ортотропію матеріалу, температурні навантаження і т. д. лише відбиває ситуацію у цій галузі зараз. Безперечно, МКЕ має суттєві переваги в порівнянні з іншими підходами та значною мірою універсальний. Разом з тим його слід розглядати як один із численних ступенів розвитку засобів чисельного дослідження при проектуванні. Загальна схема алгоритму МКЕ Метод кінцевих елементів передбачає такі основні етапи: 1. Ідеалізація фізичної системи. Під ідеалізацією розуміють процес переходу від вихідної фізичної системи до математичної моделі. Цей процес є найважливішим кроком під час вирішення технічної чи інженерної завдання. Ключовим пунктом у цьому процесі є поняття моделі, яку можна визначити як символічний пристрій, побудований для моделювання та передбачення поведінки системи. Математичне моделювання, або ідеалізація, є процесом, за допомогою якого інженер переходить від реальної фізичної системи до математичної моделі системи. Цей процес називається ідеалізацією, оскільки математичну модель необхідно абстрагувати від фізичної реальності. Як приклад реальної фізичної системи розглянемо інженерну конструкцію як плоскої пластини, навантажену поперечними силами. Математичні моделі даної системи, які інженер може використовувати для аналізу напруги в пластині, можуть бути наступними: 1) модель дуже тонкої пластини, заснована на теорії вигину мембран; 4 2) модель тонкої пластини, заснована на класичній теорії Кірхгоффа; 3) модель досить товстої пластини, заснована, наприклад, на теорії Міндліна-Рейсснера; 4) модель дуже товстої пластини, що базується на тривимірній теорії пружності. Очевидно, інженер повинен мати достатні теоретичні знання, щоб правильно вибрати відповідну математичну модель системи (конструкції), яку йому необхідно дослідити. 2. Дискретизація аналізованої області. Конструкція, що розраховується, розбивається уявними точками, лініями або поверхнями на елементи кінцевих розмірів (кінцеві елементи). Передбачається, що елементи пов'язані між собою у вузлових точках, які розташовані на їх межах. У деяких завданнях будівельної механіки як невідомі, крім вузлових переміщень, приймаються також їх приватні похідні. 3. Побудова інтерполюючих функцій. Вибирається система функцій (найчастіше шматково-поліноміальна), що однозначно визначає переміщення всередині кожного кінцевого елемента через переміщення вузлових точок. Інтерполюючі функції підбираються таким чином, щоб забезпечити безперервність шуканих величин (переміщень та їх похідних) вздовж меж елемента. 4. Виведення основних геометричних та фізичних співвідношень. На основі обраної системи інтерполюючих функцій виводяться залежності між деформаціями та переміщеннями (геометричні співвідношення), а також між напругами та деформаціями (фізичні співвідношення). 5. Побудова матриці жорсткості кінцевого елемента. За допомогою принципу Лагранжа на основі отриманих геометричних та фізичних співвідношень будується матриця жорсткості кінцевого елемента. 6. Отримання системи рівнянь методу кінцевих елементів. Кожна матриця жорсткості окремого кінцевого елемента включається до глобальної матриці жорсткості в циклі елементів. Таким чином формується система рівнянь алгебри всієї конструкції (рівняння рівноваги), яка має вигляд Kz = P , 5 де K - матриця жорсткості системи (ансамблю) кінцевих елементів; z – вектор невідомих вузлових переміщень; Р – вектор вузлових навантажень. У матриці жорсткості K, записаної вище системи рівнянь, необхідно врахувати граничні умови, оскільки інакше ця матриця буде виродженою. 7. Розв'язання системи рівнянь алгебри. Для вирішення системи лінійних рівнянь алгебри (СЛАУ) використовуються як точні, так і (при високому порядку системи) ітераційні методи. Побудовані на їх основі ефективні чисельні процедури враховують симетрію та стрічкову структуру матриці жорсткості системи. 8. Визначення деформацій та напруг. Деформації, напруги та зусилля в конструкції визначаються за допомогою знайдених вузлових переміщень на основі геометричних та фізичних співвідношень. Розглянемо деякі з цих етапів докладніше. Дискретизація аналізованої області Розбиття конструкції на кінцеві елементи - це дуже важливий крок у процедурі розв'язання задачі МКЕ, оскільки від нього багато в чому залежить точність одержуваного рішення. Успіх цьому етапі забезпечують, насамперед, наявні інженерні навички. Невдало виконане розбиття області на кінцеві елементи може призвести до помилкових результатів. При призначенні сітки КЕ постає завдання оптимального розбиття конструкції на підобласті. При цьому слід мати на увазі, що розміри елементів повинні бути достатньо малими для того, щоб забезпечити прийнятну точність рішення, з іншого боку, використання густої сітки призводить до великих систем рівнянь алгебри, рішення яких пов'язане зі значним обсягом обчислювальної роботи. У процесі розбиття області на кінцеві елементи необхідно враховувати деякі загальні уявлення про остаточні результати розрахунку для того, щоб зменшити розміри кінцевих елементів у зонах концентрації напруги, де шукані величини швидко змінюються, і збільшити розміри КЕ там, де шукані величини змінюються повільно. Важливим моментом у процесі розв'язання задачі МКЕ є нумерація вузлів сітки, оскільки від цього залежить ширина стрічки 6 матриці роздільних рівнянь, відповідно час рахунку та обсяг використовуваної пам'яті ЕОМ. В даний час розроблено сервісні програми автоматизованого розбиття конструкції на кінцеві елементи та раціональну нумерацію вузлів. Побудова інтерполюючих функцій МКЕ заснований на апроксимації безперервної функції, визначеної по всій області, дискретною моделлю за допомогою шматково-безперервних функцій, визначених на підобластях (кінцевих елементах). Запишемо переміщення, що є функціями координат довільної точки кінцевого елемента через компоненти вектора вузлових переміщень за допомогою інтерполюючої функції (функції форми або базисної функції): u = Nz , (1) де N = [ N1 N 2 … N s ] - матриця функції форми ; z = ( z1 z2 ... zs ) - Вектор вузлових переміщень кінцевого елемента (КЕ); s – кількість ступенів свободи КЕ. Функції (1) повинні відповідати критеріям повноти та спільності. Розглянемо їх. 1. Критерій повноти. Інтерполіруюча функція повинна забезпечувати постійні значення величин, що розглядаються при зменшенні розмірів елемента. Для виконання цієї умови інтерполююча функція повинна бути повним поліном як мінімум ступеня р, де р - найвищий порядок похідної, що входить у функціонал. T Умова повноти задовольняється у тому випадку, коли s ∑ Ni = 1 . i = 1 2. Критерій спільності. Інтерполююча функція повинна бути безперервна разом зі своїми похідними до (n – 1)-го порядку включно (де n - максимальний порядок похідної в підінтегральному вираженні функціоналу енергії) на межі між елементами. Критерії повноти та спільності є достатні умови збіжності методу кінцевих елементів. При їх виконанні зі зменшенням розмірів кінцевого елемента наближені рішення МКЕ 7 монотонно сходяться до точного рішення. Сказане зовсім на означає, що порушення цих критеріїв призводить до неможливості отримання достовірного рішення. Існують несумісні та навіть неповні елементи, які дають високу точність та швидку збіжність. Висновок основних геометричних та фізичних співвідношень У загальному вигляді залежність між деформаціями та переміщеннями (геометричні співвідношення) записується так: ε = Bz , (2) де ε - вектор деформацій; z – вектор вузлових переміщень; - матриця, що зв'язує вектор вузлових переміщень з вектором, що містить компоненти тензора деформації. Таким чином, вважається, що залежність (2) між деформаціями та вузловими переміщеннями має лінійний характер. Лінійній залежності відповідають такі умови роботи конструкції, коли деформації та кути повороту малі порівняно з одиницею, а квадрати кутів повороту малі порівняно з відповідними компонентами деформації. Фізичні співвідношення, що визначають залежність між напругами та деформаціями, мають вигляд σ = Dε, (3) де σ - вектор, що містить компоненти тензора напруги; D – матриця пружності. Рівняння стану (3) є узагальненим законом Гука, що встановлює прямо пропорційну залежність між напругами і деформаціями, справедливий для певного класу матеріалів на певній ділянці графіка σ − ε . Побудова матриці жорсткості кінцевого елемента Розв'язання задач розрахунку конструкцій базується на двох основних підходах. У першому випадку вирішуються диференціальні рівняння із заданими граничними умовами. У другому випадку записується умова стаціонарності інтегральної величини, пов'язаної з роботою напруг і зовнішнього прикладеного навантаження і являє собою повну потенційну енергію системи. Для розрахунку конструкцій у межах МКЕ використовується другий підхід. Як відомо, повна потенційна енергія пружної системи визначається за формулою 8 (z) = W(z) − A(z) , де W - потенційна енергія деформації; А – потенціал зовнішніх сил. Потенційна енергія деформації пружної системи визначається співвідношенням W= 1 T ε σ dV , 2V де V - об'єм, який займає тіло, а потенціал зовнішніх розподілених навантажень визначається за формулою A = uT p dS , S де р - вектор зовнішніх розподілених навантажень; S – площа, на якій прикладено навантаження. При цьому деформації та напруги, що входять до формули для потенційної енергії, виражаються через вузлові переміщення. Отримання рівнянь МКЕ в переміщеннях засноване на одному з фундаментальних енергетичних принципів механіки - принципі Лагранжа, згідно з яким для системи, яка перебуває в стані рівноваги, повна потенційна енергія набуває стаціонарного значення. Ця умова записується у вигляді ∂Π = 0. ∂z Будемо вважати, що значення повної потенційної енергії для всієї області V дорівнює сумі енергій окремих кінцевих елементів: m () m (() ()) Π (z) = ∑ Π z = ∑ W i z i − Ai z i , i =1 i i i =1 (4) де m - кількість кінцевих елементів. Тоді ∂Π m ⎛ ∂W i (z) ∂Ai (z) ⎞ = ∑⎜ − ⎟ = 0. ∂z i =1 ⎝ ∂z ∂z ⎠ (5) Розглянемо окремий кінцевий елемент, індекс i при цьому 1 T 1 T (Bz)T DBz dV − ∫ (Nz)T p dS = ε σ dV − u p dS = ∫ ∫ ∫ 2V 2V S S 1 1 = zT (∫ BT DBz dV) z − zT ∫ N T p dS zT Kz − zTP, (6) 2 2 S Π(z) = де K = ∫ BT DB dV - матриця жорсткості кінцевого елемента, а (7) V TP = ∫ N p dS - вектор вузлових навантажень. (8) S Отримання системи рівнянь методу кінцевих елементів Для виконання операції підсумовування необхідно перетворити вектори вузлових переміщень z(i) та вузлового навантаження P(i) для окремого кінцевого елемента у відповідні вектори z та Р для всієї системи, що може бути зроблено за допомогою деякої булевої матриці H(i), що містить як елементи тільки нулі і одиниці: z (i) = H (i) z; P(i) = H(i) P. (9) Підстановка формул (9) у вираз для повної потенційної енергії кінцевого елемента (6) дає: () () T 1 (i) T (i) (i) H z K H z − H (i) z H (i ) P = 2 T T 1 = zT H (i) K (i) H (i) z − zT H (i) H (i) P. 2 Тоді диференціювання по z, відповідно до формули (5), призводить до системи рівнянь: Π (i) = m ∑ i = 1 (T T) H (i) K (i) H (i) z − H (i) H (i) P = 0 , (10) де використано правило диференціювання матричних співвідношень ∂ T z Kz = 2 Kz. ∂z Система (10), яку можна записати у вигляді () Kz = P , (11) 10 являє собою систему лінійних рівнянь алгебри методу кінцевих елементів, що є рівняннями рівноваги в переміщеннях. Як правило, рішення системи (11) виконується методом Гауса чи ітераційними методами. (i)T Матриця жорсткості окремого елемента H K (i) H (i) , що фігурує у формулі (10), являє собою розширену матрицю, розмірність якої дорівнює розмірності глобальної матриці. Тому використання процедури підсумовування у формулі (10) при чисельній реалізації МКЕ є неефективним. У практичних розрахунках виконується пряме побудова глобальної матриці жорсткості. У цьому випадку будується матриця для окремого кінцевого елемента за формулою (7), що має розмірність S×S. Потім рядкам і стовпцям цієї матриці приписуються номери глобальних ступенів свободи, що дозволяє визначити розташування коефіцієнтів матриці жорсткості КЕ в глобальній матриці жорсткості. Після цього попередньо обнулену глобальну матрицю заносяться коефіцієнти матриці жорсткості КЕ те місце, що визначено їх адресою. Нехай, наприклад, система складається з двох КЕ, що містять по два вузли, у кожному з яких є одне невідоме (один ступінь свободи). Загальна кількість вузлів системи – 3, розмірність глобальної матриці 3×3, елементи з'єднані між собою у 2-му вузлі. Матриці жорсткості 1-го і 2-го КЕ, з відповідною нумерацією коефіцієнтів, і глобальна матриця мають вигляд K(1) = 1 k11 1 1 k21 1 k12 ⎥; 1 k22 ⎦⎥ K (2) = 2 ⎡ k22 ⎢ 2 ⎢⎣ k32 1 1 ⎡ k11 k12 2 ⎤ ⎢ 1 k23 1 2 = ; K ⎥ ⎢ k21 k22 + k22 2 k33 ⎥⎦ ⎢ 2 k32 ⎢⎣ 0 0 ⎤ ⎥ 2 k23 ⎥. 2 ⎥ k33 ⎥⎦ У матриці K системи рівнянь (11) необхідно врахувати граничні умови, інакше вона буде виродженою, тобто її визначник дорівнюватиме нулю. Облік граничних умов може бути здійснено трьома різними способами 1. З матриці K видаляються k рядок і k стовпець, відповідні переміщенню zk = 0 . Після цього рядки та стовпці матриці 11 перенумеровуються. Відповідно, зменшується розмірність вектора вузлових переміщень. 2. Рівняння zk = 0 , що відповідає граничній умові, формується у складі матриці K. Для отримання zk = 0 у матриці K k-й рядок і k-й стовпець, а також відповідний елемент у векторі зовнішніх навантажень P заповнюються нулями. На місце діагонального елемента rrr у матриці K ставиться одиниця. В результаті порядок матриці не змінюється, а задані переміщення набувають нульових значень. 3. Для отримання zk = 0 діагональний елемент rrr множиться на велике число. Порядок матриці у своїй не змінюється. Визначення деформацій та напруг У результаті розв'язання системи рівнянь (11) визначається вектор вузлових переміщень усієї конструкції. На основі знайдених значень вузлових переміщень за формулою (2) визначається вектор деформацій КЕ, а за формулою (3) вектор напруг. Двовимірний симплекс-елемент Класифікація кінцевих елементів може бути проведена відповідно до порядку багаточленів – функцій цих елементів. При цьому розглядаються три групи елементів: 1) симплекс-елементи; 2) комплекс-елементи; 3) мультиплекс-елементи. Симплекс-елементам відповідають багаточлени першого ступеня. Комплекс-елементам – багаточлени вищого порядку. У симплекс-елементі число вузлів дорівнює розмірності простору +1. У комплекс-елементі число вузлів більше за цю величину. Для мультиплекс-елементів також використовуються багаточлени високого порядку, але межі елементів повинні бути паралельні координатним осям. Розглянемо формування матриці жорсткості для двовимірного симплекс-елемента. Двовимірний симплекс-елемент є трикутником з вузлами, розташованими в його вершинах (рис. 1). 12 Мал. 1. Двовимірний симплекс-елемент Вузли КЕ нумеруються проти годинникової стрілки, починаючи з деякого довільно обраного i-го вузла. Координати i-го, j-го і k-го вузлів по осі x позначені через xi, x j, xk, по осі y - через yi, y j, yk. Кожен вузол має два ступені свободи - переміщення u вздовж осі x і переміщення v вздовж осі y. Інтерполюючі функції, що визначають переміщення довільної точки КЕ вздовж осей х і у, приймаються у вигляді u(x, y) = α 0 + α1 x + α 2 y; (12) v(x, y) = α3 + α 4 x + α5 y. Коефіцієнти α 0 ,…, α5 визначаються за допомогою граничних умов: u = ui , v = vi при x = xi та y = yi ; u = u j , v = v j при x = x j і y = y j; u = uk, v = vk при x = xk та y = yk. Визначимо коефіцієнти α0, α1, α2. Для цього підставимо граничні умови для функції і в перший вираз (12), що призведе до системи рівнянь: ui = α 0 + α1 xi + α 2 yi; u j = α 0 + α1 x j + α 2 y j; uk = α 0 + α1 xk + α 2 yk. 13 ⎡1 xi ⎢ Або ⎢1 x j ⎢1 x k ⎣ yi ⎤ ⎧α 0 ⎫ ⎧ ui ⎫ ⎥⎪ ⎪ ⎪ ⎪ y j ⎥ ⎨ α1 . ⎪ ⎪ ⎪ ⎪ yk ⎥⎦ ⎩α 2 ⎭ ⎩uk ⎭ (13) Визначник системи (13) дорівнює подвоєній площі F трикутного елемента: 1 xi 1 xj yi y j = 2F . 1 xk yk (14) Тоді за правилом Крамера α0 = ui uj xi xj yi yj uk xk yk 1 xi 1 xj yi yj 1 xk yk = ui uj xi xj yi yuk xk yk 2F або 1 ⎡ ui x k yk u j (xi yk − x k yi) − uk xi y j − x j yi ⎤ . ⎣ ⎦ 2F Аналогічно визначаються коефіцієнти α1 та α2. Після підстановки виразів для α0 , α1 , α 2 у першу формулу (12) маємо 1 ⎡ u(x, y) = (ai + bi x + c i y) ui + a j + bj x + c j yu j + 2F ⎣ + (ak + bk x + c k y) uk ⎤⎦ , (15) α0 = () (()) де ai = x j yk − x k y j ; bi = y j − yk; ci = xk − x j. (16) Інші коефіцієнти (15) виходять за формулами (16) циклічною перестановкою індексів (індекс i замінюється на індекс j, індекс j - на індекс k, індекс k - на індекс i). Аналогічно: v(x, y) = 1 ⎡ (ai + bi x + c i y) vi + a j + bj x + c j y v j + 2F ⎣ + (ak + bk x + c k y) vk ⎤⎦ . () 14 (17) Тоді в матричному вигляді ⎡ Ni ⎧u ⎫ u = ⎨ ⎬ = Nz = ⎢ ⎩v ⎭ ⎢⎣ 0 0 Nj 0 Nk Ni 0 Nj 0 ⎧ ui ⎫ ⎪ ⎪ u j ⎪⎪ ⎥ ⎨ ⎬, N k ⎥⎦ ⎪ v j ⎪ ⎪u ⎪ ⎪ k⎪ ⎩⎪ vk ⎭⎪ 1 1 a j + bj x + c j y ; (ai + bi x + c i y); N j = 2F 2F (18) 1 Nk = (ak + bk x + c k y). 2F Геометричні співвідношення, що пов'язують деформації та переміщення в рамках плоскої задачі теорії пружності, записуються з використанням формул (15), (17) наступним чином: (де Ni =) ∂u 1 ∂v 1 bi ui + b j u j + bk uk ; ε y = ci vi + c j v j + ck vk; = = ∂x 2 F ∂y 2 F ∂u ∂v 1 ci ui + c j u j + ck uk + bi vi + b j v j + bk vk γ xy = + = ∂y ∂x 2 F або у матричному вигляді: εx = () ((⎧ε ⎫ ⎡bi x ⎪ ⎪ 1 ⎢ ε = ⎨ εy ⎬ = ⎢0 ⎪ ⎪ 2F ⎢ ⎩ γ xy ⎭ ⎣ci ⎡bi 1 ⎢ де B = ⎢ ) ci 0 cj 0 bi cj bj ck 0 bj 0 bk ci 0 cj 0 bi cj bj ck ⎧ ui ⎫ ⎪v ⎪ 0 ⎤⎪ i ⎪  ⎪ bk ⎦ ⎪ ⎪ u ⎪ k⎪ ⎩⎪ vk ⎭⎪ 0⎤ ⎥ ck ⎥ – матриця градієнтів. ⎥ bk ⎦ 15 (19) Величина подвоєної площі кінцевого елемента 2F у виразі (19) підраховується за формулою (14). Фізичні співвідношення, що визначають залежність між напругами і деформаціями в плоскій задачі теорії пружності, записуються у вигляді σ = Dε , ⎡ ⎤ ⎢ 1 ν1 0 ⎥ ⎥ E1 ⎢ де D = 1 0 ν 1 ⎥ - матриця упр 2 ⎢ 1 − ν1 ⎢ 1 − ν1 ⎥ ⎢0 0 ⎥ 2 ⎦ ⎣ (20) У разі плоскої деформації (ε z = 0) у формулі (20) слід прийняти E1 = E ν ; ν = , 1 1− ν 1 − ν12 (21) а для узагальненого плоского напруженого стану (σ z = 0) E1 = E ; ν1 = ν. (22) Формули (21) та (22) відповідають ізотропному матеріалу з модулем пружності Е та коефіцієнтом Пуассона ν . Не представляє великої праці побудувати матриці пружності для ортотропного матеріалу, коли жорсткі характеристики різні у двох взаємно перпендикулярних напрямках. Оскільки матриці і D містять тільки константи, об'ємний інтеграл, що визначає матрицю жорсткості елемента у формулі (7), легко обчислюється: K = BT DB dV = BT DB dV V (23) V або K = BT DB tF . (24) У формулі (24) t - товщина елемента; F – площа елемента. Зазвичай матриця жорсткості (24) визначається чисельно. Для цього спочатку знаходяться чисельні значення коефіцієнтів матриць В і D, а потім виконується перемноження відповідно до виразу (23) або (24). ПРИКЛАД ВИКОНАННЯ ЛАБОРАТОРНОЇ РОБОТИ Перед виконанням лабораторної роботи ще раз нагадаємо основні етапи методу МКЕ: 1. Пружне тіло розбивається на елементи. Об'ємне тіло на тетраедри або паралелепіпеди. Плоске тіло - на трикутники та прямокутники. 2. До кожного елемента складається матриця жорсткості з допомогою функції форми. Функція форми є способом апроксимації невідомої функції переміщень. 3. Матриці жорсткості елементів поєднуються в єдину матрицю жорсткості для всього тіла. 4. Вирішуючи систему рівнянь, знаходять вузлові переміщення. 5. За допомогою рівнянь теорії пружності визначаються деформації та напруження у вузлових точках тіла. У наведеному прикладі вирішується плоска задача теорії пружності. Кільце, навантажене двома силами (рис. 2, в) має дві осі симетрії, тому для підвищення точності розрахунку розглянемо одну четверту частину кільця (рис. 3). На осях симетрії мають виконуватися граничні умови рівності нулю переміщень, перпендикулярних до осей симетрії. Розглянуту чверть кільця розбиваємо на кінцеві трикутні елементи (див. рис. 2, б). Трикутний елемент має 6 ступенів свободи (незалежних вузлових переміщень). Нумерація вузлових переміщень елементі починається в нижньому лівому вузлі трикутника і продовжується проти годинникової стрілки. Горизонтальні переміщення – непарні, вертикальні – парні. Нумерація вузлів всього тіла та кінцевих елементів – по стовпцях зверху вниз, зліва направо. Розміри елементів можуть бути різними (що менше елемент, то вища точність розрахунків). У нашому прикладі всього 66 вузлів та 100 кінцевих елементів. Положення розрахованих вузлів показано на рис. 3, б. Розрахунок координат вузлів наведено на рис. 4. Введення координат вузлів - це копітка праця, і при великій кількості вузлів краще автоматизувати цю роботу. 17 а б Рис. 2. Схема навантаження і трикутний кінцевий елемент Рис. 3. Розрахункова схема та координати вузлів На рис. 4 наводиться розрахунок полярних координат вузлів та їх перетворення на прямокутні (декартові) координати. Тут r1 і r2 - зовнішній та внутрішній радіуси кільця; t – товщина кільця; φ1 і φ2 - початкове та кінцеве значення кутової координати; X0 та Y0 - декартові координати полюса (початки полярних координат); nr і nφ - число вузлів у стовпці (вздовж радіусу) і в ряду (по куту охоплення частини тіла, що розглядається). Результати розрахунків координат вузлів наведено на графіку (рис. 5). 18 Мал. 4. Розрахунок координат вузлів елементів Мал. 5. Сітка вузлів 19 Кропіткою та трудомісткою є також завдання складання матриці індексів. На рис. 6 наведена програма складання матриці індексів, що використовується в прикладі. Тут же наведено автоматичний розрахунок граничних умов і, залежно від кількості вузлів, номерів переміщень, у яких переміщення осях симетрії дорівнює нулю. Автоматизація розрахунків координат вузлів, матриці індексів та граничних умов дозволяє для даної схеми змінювати кількість вузлів (рис. 7). Мал. 6. Програма розрахунку матриці індексів Рис. 7. Матриця індексів, координати вузлів, номери та задані переміщення Координати вузлів, матрицю індексів та граничні умови можна записати в окремі файли, як це показано на рис. 8. Надалі ці файли можна вважати за допомогою функції READPRN та використовувати в іншому документі. 21 Мал. 8. Запис координат, матриці індексів та граничних умов у зовнішні файли Даний розрахунок проведено з урахуванням розмірності (рис. 4). Облік розмірностей вносить додаткові труднощі у досить складний розрахунок, особливо у введенні матриць від розмірних величин (рис. 9). Мал. 9. Вектор сил і матриця індексів переміщень 22 кожного трикутного елемента 3 вузла і 6 вузлових переміщень. Матриця індексів переміщень (матриця зв'язку глобальних номерів вузлових переміщень тіла з локальними номерами вузлових переміщень елементів) отримана шляхом подвоєння MIU матриці. Матриця жорсткості елемента розраховується за формулою K = BT DB dV = BT DB dV . V (23) V Тут B = ∂T N , де D - матриця внутрішньої жорсткості, що містить пружні постійні матеріалу E , ; ∂ - матричний диференціальний оператор, що означає певну послідовність присвоєння знака диференціювання; N – матриця функцій форми. Для трикутного елемента функція форми - рівняння площини, що визначається виразами (18). Як було показано вище (див. вираз (19)), матриця B = ∂T N містить константи, які залежать лише від координат вузлів. Наведемо розрахунок коефіцієнтів на формування матриці жорсткості елемента (рис. 10) і розрахунок площі елементів (рис. 11). Мал. 10. Розрахунок коефіцієнтів на формування матриці жорсткості елемента 23 Рис. 11. Розрахунок площі елементів Матриця внутрішньої жорсткості D наведена нижче (рис. 12) і записана як умовного оператора: різні матриці для плоского напруженого NDS = 0 і плоского деформованого стану NDS = 1. Рис. 12. Формування матриці внутрішньої жорсткості Для трикутного елемента інтеграл за обсягом дорівнює добутку підінтегрального виразу обсяг. Формула (23) для розрахунку матриці жорсткості елемента наведена зверху. Матриця жорсткості системи формується з допомогою матриці індексів. 24 Облік граничних умов супроводжується перебудовою матриці жорсткості системи та вектора сил (рис. 13). Мал. 13. Формування матриці жорсткості системи та облік граничних умов Вузлові переміщення визначаються шляхом обігу матриці жорсткості. Мал. 14. Визначення переміщень вузлів, деформацій та напруг у центрі елемента Відповідно до рівнянь теорії пружності ε = ∂T u , де u - вектор переміщень. Відповідно до рівнянь зв'язку вузлових переміщень ∆ та переміщень u довільної точки u = N ∆ . 25 () Звідси деформація елемента ε = ∂T N ∆ . З фізичних рівнянь теорії пружності (закон Гука) напруги = Dε. Складність розрахунку полягає в акуратному використанні індексів елементів, вузлів, стовпців, рядків, присвоєння індексів значень, взятих із матриці індексів. Для трикутного елемента функція форми лінійна, тому похідні функції форми, деформації і напруги, знайдені на рис. 14, постійні у всій площі елемента. Напруги у вузлах тіла визначаються як середня арифметична напруга або деформація у всіх елементах, що сходяться у вузлі. Розрахунок напруг та деформацій у вузлах тіла наведено на рис. 15. Мал. 15. Визначення переміщень вузлів тіла, напружень та деформацій у центрі кожного елемента 26 На рис. 15 показано визначення 4-го значення деформації та напруги в кожному елементі, що не враховуються в матриці внутрішньої жорсткості D. Працюючи з документом Mathcad, якщо ми опустимо вираз NDS: = 1 нижче виразу NDS: = 0, то ми можемо побачити результати розрахунку вже не при плоскому напруженому стані, а за плоскої деформації. Результати розрахунку наведено на рис. 16. Мал. 16. Результати расчета 27 Завдання до лабораторної работе Вирішити методом кінцевих елементів (МКЭ) плоску завдання теорії пружності. Кільце, навантажене двома силами (рис. 3, б), має дві осі симетрії, тому підвищення точності розрахунку необхідно розглянути одну четверту частину кільця. За допомогою відповідної сітки вона поділяється на систему кінцевих трикутних елементів. Число вузлів уздовж радіусу nr і число вузлів по куту охоплення частини тіла nφ, що розглядається, вибираються з таблиці індивідуальних завдань відповідно до варіанту. Визначити деформації та напруження при плоскому напруженому стані та при плоскому деформованому стані. Варіанти індивідуальних завдань № Число вузлів варі- вздовж радіуанта са, nr 1 8 2 9 3 10 4 11 5 12 6 13 7 8 8 9 9 10 10 11 11 12 12 13 13 8 14 9 15 4 5 6 7 8 9 5 6 7 8 9 4 6 7 8 № варианта 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 Число узлов вдоль радиуса, nr 11 12 13 9 10 11 12 13 10 11 12 13 11 12 13 Число вузлів по куту охоплення, nφ 9 4 5 6 4 8 9 7 5 6 7 6 4 5 6 Зміст звіту У робочому каталозі студента мають бути створені два файли, що містять налагоджені документи в системі Mathcad, що відповідають розрахункам при плоскості та при плоскому деформованому стані. 28 Звіт з лабораторної роботи повинен містити: 1) назву лабораторної роботи; 2) мету лабораторної роботи; 3) завдання; 4) скопійовані з екрана монітора налагоджені документи виконання завдання. Контрольні питання 1. Навіщо призначений метод кінцевих елементів? 2. Навести загальну схему розрахунку методом кінцевих елементів. 3. Що таке ідеалізована конструкція? 4. З якої системи лінійних рівнянь алгебри визначаються переміщення вузлів? 5. Як визначається вид апроксимуючого полінома? 6. Які функції називаються базисними функціями? 7. Як визначаються функції форми? 8. Як визначаються функції переміщень? 9. Як визначаються функції деформацій? 10. Як визначається напруга? 11. Який принцип використовується визначення узагальнених сил? 12. Як визначається матриця жорсткості елемента? 13. Як складається матриця жорсткості системи? 14. Від урахування яких впливів утворюються вільні члени системи канонічних рівнянь? 15. Чим можуть бути спричинені попередні деформації? 16. Чим можуть бути викликані попередні напруги? 17. Як визначаються реактивні зусилля внаслідок окремих впливів? 18. Як у системі Mathcad здійснюється розрахунок координат вузлів сітки? 29 Бібліографічний перелік 1. Дарков В.А. Будівельна механіка: навч. для будує. спец. вузів/В.А. Дарків, Н. Н. Шапошников. Вид. 8-е, перероб. та дод. М.: Вищ. шк., 1986. 607 с. 2. Макаров Є.Г. Інженерні розрахунки в Mathcad 14: навчальний курс. СПб. : Пітер, 2007. 592 с. 3. Трушін С.І. Метод кінцевих елементів. Теорія та завдання: навчальний посібник. М.: Вид-во АСВ, 2008. 256 с. 4. Хечумов Р.А. Застосування методу кінцевих елементів до розрахунку конструкцій/Р.А. Хечумов, Х. Кеплер, В.І. Прокіпів. М.: Вид-во АСВ, 1994. 353 с. 30 Навчальне видання Ігнатьєв Олександр Володимирович, Михайлова Наталія Анатоліївна, Єрещенко Тетяна Володимирівна МЕТОД КІНЦЕВИХ ЕЛЕМЕНТІВ ТА ЙОГО РЕАЛІЗАЦІЯ У СЕРЕДОВИЩІ MATHCAD Лабораторний практикум з дисципліни «Аналітичні та чисельні методи вирішення рівнів. РІО О.Є. Горячева Зав. редакцією М.Л. Піщана Редактор О.А. Шипунова Комп'ютерна правка та верстка Н.А. Деріна Підписано до друку 30.06.10. Формат 60х84/16. Папір офсетний. Друк трафаретний. Гарнітура Таймс. Ум. піч. л. 1,9. Уч.-вид. л. 1,7. Тираж 100 екз. Замовлення №70 Державна освітня установа вищої професійної освіти «Волгоградський державний архітектурно-будівельний університет» Редакційно-видавничий відділ Сектор оперативної поліграфії ЦІТ 400074, Волгоград, вул. Академічна, 1 31