Чисельні методи у механіці

Вступ. У даному розділі будуть розглянуті особливості формування системи лінійних алгебраїчних рівнянь методу скінченних елементів (СЛАР МСЕ) на прикладі тривузлового трикутного елемента з лінійною інтерполяцією переміщень, який використовується для розв'язання плоскої задачі теорії пружності. Скорочено називатимемо такий елемент лінійним трикутним скінченним елементом.

Цей елемент має ряд відмітних особливостей:

1) він належить до сім'ї так званих ізопараметричних елементів, про що йдеться в наступних розділах;

2) він дозволяє одержати вирази елементних матриць жорсткості і елементних векторів сил у замкненій формі, що означає відсутність необхідності в чисельному інтегруванні при обчисленні елементних матриць жорсткості і елементних векторів сил;

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

На завершення хотілося б відзначити, що лінійний трикутний кінцевий елемент має певне історичне значення. Він був одним із двох перших скінченних елементів, які були представлені у статті Мартіна, Тернера, Клоха і Топпа в 1956 році. Ця публікація загальновизнано вважається початком сучасного методу скінченних елементів.

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

Параметричне подання функцій. У подальшому викладенні істотним чином використовуватиметься поняття параметричного подання функції. Разом з алгоритмами чисельного інтегрування даний підхід виконує ключову роль у системному розробленні різних типів скінченних елементів для розв'язання двовимірних і тривимірних задач механіки деформованого тіла.

Основна ідея параметричного подання функції може бути продемонстрована на простому прикладі. Розглянемо канонічне рівняння кола одиничного радіуса з центром на початку координат:

. (3.53)

Окрім канонічного подання (3.53), рівняння кола може бути подане у вигляді стандартної функціональної залежності

, (3.54)

а також в полярній системі координат:

,

. (3.55)

У поданні (3.55) є полярним кутом, який може бути інтерпретований у вигляді незалежного скалярного параметра, який є змінним у заданих межах. Очевидно, що підстановка (3.55) в (3.53) дає тотожну рівність.

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

,

, (3.56а)

,

,

, (3.56б)

де - незалежний скалярний параметр, який змінюється у заданих межах.

Вигляд функцій, що входять до (3.56), може бути будь-яким. Головна умова, яка накладена на ці функції, полягає у тому, щоб кожному значенню параметра однозначно відповідала точка на площині, яка має координати , або в просторі - .

Аналогічно визначається параметричне подання шматка площини на координатній площині або шматка поверхні у тривимірному просторі:

,

, (3.57а)

,

,

. (3.57б)

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

Параметричне подання функцій істотно використовується в алгоритмах методу скінченних елементів.

Система координат трикутного елемента. Геометрія тривузлового трикутного скінченного елемента задається трьома кутовими точками в площині - вузлами скінченного елемента (рис. 3.13). Локальні номери вузлів на елементі - 1, 2, 3, причому вузли нумеруються так, щоб обхід контуру елемента був проти годинникової стрілки, якщо дивитися з боку зовнішньої нормалі до площини пластини. Нагадаємо, що зовнішня нормаль до площини пластини збігається з додатним напрямом глобальної осі z. Положення вузлів визначається декартовими координатами у глобальній системі координат (3.58а).

Рисунок 3.13 - Лінійний трикутний скінченний елемент:

(а) - геометрія елемента; (б) - додатний напрямок обходу контура елемента

Лінійний трикутний скінченний елемент має шість ступенів вільності, які утворені компонентами елементного вектора вузлових переміщень (3.58б).

Позначимо площу елемента через А. Згідно з відомою формулою лінійної алгебри площа плоского трикутника може бути обчислена за допомогою визначника, який складений з координат вершин трикутника. Отже, площа лінійного трикутного скінченного елемента дорівнює

. (3.59)

Зазначимо, що площа області, яка визначається за формулою (3.59), додатна, якщо локальна нумерація вузлів (1, 2, 3) виконана проти годинникової стрілки, якщо дивитися з боку зовнішньої нормалі до площини пластини, як це буде відмічено вище і показано на рисунку 3.13. Інакше знак А буде від'ємним, що свідчить про помилку. Така угода використовуватиметься і в подальшому викладенні матеріалу.

Окрім глобальної системи координат , введемо локальну параметричну систему координат у площині елемента, яка задається трьома змінними величинами . (3.60)

Ці три локальні параметричні координати в літературі з МСЕ називаються також природними і трикутними. За допомогою введеної системи координат встановлюється взаємооднозначна відповідність між глобальними декартовими координатами довільної точки елемента і трьома скалярними параметрами :

(3.61)

Відразу звернемо увагу на уявну суперечливість між формулами (3.61) і (3.57а), яка легко усувається, якщо допустити наявність зв'язку між параметрами . Насправді, так воно і є: три трикутні координати не є незалежними, а зв'язані одним співвідношенням

. (3.62)

Очевидно, що із співвідношення (3.62) можна виразити, наприклад, координату через інші дві . Тоді співвідношення (3.61) набуде стандартного вигляду (3.57а). Проте використання трьох параметрів виявляється зручнішим з погляду обчислень і розроблення ефективних чисельних алгоритмів, про що завжди необхідно пам'ятати, коли маєш справу з обчислювальною механікою.

На рисунку 3.14 дано геометричне пояснення математичного значення введених координат.

Рисунок 3.14 - Трикутні координати

Три рівняння, записані у вигляді , (3.63)

геометрично є множиною прямих ліній, паралельних стороні, протилежній i-му вузлу.

Наприклад, рівняння сторін трикутного елемента 1-2, 2-3 і 3-1 матимуть такий вигляд:

сторона 1-2: , сторона 2-3:, сторона 3-1:. (3.64а)

Три вершини трикутного елемента матимуть такі локальні координати:

вузол 1: , вузол 2: , вузол 3: . (3.64б)

Як приклад, наведемо також координати середніх точок сторін трикутного елемента і геометричного центра елемента:

середній вузол сторони 1-2: , середній вузол сторони 2-3:, середній вузол сторони 3-1:, центр мас елемента:. (3.64в)

Ще раз нагадаємо, що три параметри не є незалежними і зв'язані одним додатковим співвідношенням

.

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

, (3.65)

де - довільні постійні коефіцієнти.

Нехай тепер ми бажаємо ввести апроксимацію даної функції в межах трикутного елемента. Позначимо її через , (3.66).

значення даної функції у відповідних вузлах скінченного елемента - так звані вузлові значення функції.

За допомогою цих значень можуть бути визначені невідомі коефіцієнти у виразі (3.65). Насправді, одержимо:

(3.67)

Співвідношення (3.67) являють собою систему трьох лінійних алгебраїчних рівнянь щодо невідомих коефіцієнтів . Розв'язавши систему (3.67), одержимо однозначний вираз для функції (3.65) через її вузлові значення. Проте це вираз достатньо складний і незручний у використанні. Тому використовують подання функції не в глобальній, а в локальній системі координат на елементі. Очевидно, що якщо встановлена взаємооднозначна відповідність між глобальною і локальною системами координат, то і вирази функції в різних координатах будуть відповідними один до одного.

Використовуючи даний підхід, подамо функцію (3.65) у трикутних координатах через її вузлові значення таким чином:

(3.68)

Співвідношення (3.68) називається лінійним співвідношенням інтерполювання довільної безперервної функції f у трикутних координатах.

Перетворення координат. Математичні величини, які тісно пов'язані з геометрією трикутного скінченного елемента, природним чином виражаються через трикутні координати . З іншого боку, такі фізичні величини, як переміщення, деформації і напруження часто можуть бути виражені в декартових координатах . Отже, необхідні рівняння, які зв'язують ці дві системи координат і за допомогою яких можна переходити від однієї системи координат до іншої.

Використовуючи універсальне співвідношення (3.68), глобальні декартові координати довільної точки на трикутному елементі можуть бути виражені через глобальні координати вузлів елемента і локальні трикутні координати цієї самої точки :

, (3.69)

де - глобальні декартові координати вузлів трикутного елемента.

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

,

. (3.70)

Розв'яжемо систему (3.69) стосовно вектора . Знайдемо зворотну матрицю до матриці системи (3.69) і запишемо розв'язання у вигляді

, (3.71)

де А - визначник матриці системи (3.69), який чисельно збігається з площею трикутного елемента.

Введемо позначення:

(3.72)

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

.

У позначеннях (3.72) розв'язання (3.71) набуде вигляду

. (3.73)

Обчислення частинних похідних. Із співвідношень (3.70) і (3.73) легко обчислюються частинні похідні за відповідними координатами:

(3.74)

У другому рядку формул (3.74) послідовність є круговою перестановкою індексів. Наприклад, якщо , то і . Якщо , то і .

Розглянемо довільну безперервну функцію двох змінних, яка змінюється лінійно в межах деякої плоскої трикутної області. Ця функція може бути представлена як в глобальній декартовій системі координат , так і в локальній трикутній системі координат . Представимо її як складну функцію багатьох змінних:

. (3.75)

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

(3.76)

Підставимо (3.74) в (3.76). Одержимо вирази для часткових похідних довільної безперервної функції за глобальними координатами через три часткові похідні тієї самої функції за локальними трикутними координатами:

,

. (3.77)

Остаточно в матричній формі система часткових похідних довільної безперервної функції за глобальними координатами набуде вигляду

. (3.78)

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

Рисунок 3.15 - Інтерполяція переміщень в межах трикутного лінійного скінченного елемента

3.8 Виведення рівнянь трикутного лінійного скінченного елемента

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

(3.79)

Тут - компоненти вектора переміщень у вузлах скінченного елемента (рис. 3.15).

Інтерполяційне співвідношення для переміщень (3.79) може бути також подане в стандартній для методу скінченних елементів матричній формі:

, (3.79а)

де - матриця функцій інтерполювання; - елементний вектор вузлових переміщень.

Очевидно, що в даному випадку трикутного лінійного скінченного елемента функції інтерполяції дорівнюють трикутним координатам

. (3.80)

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

, (3.81)

де - матриця градієнтів.

Обчислимо матрицю градієнтів

(3.82)

Функції інтерполяції можуть бути інтерпретовані як складні функції багатьох змінних

(3.83)

З іншого боку, в явному вигляді вони збігаються з трикутними координатами елемента (див. формулу 3.80). За допомогою формул для обчислення часткових похідних функцій подібного типу за глобальними координатами (3.77) і (3.74) одержимо вираз матриці градієнтів

. (3.84)

Звернемо увагу, що деформації (3.81) вийшли постійними на елементі, оскільки компоненти матриці градієнтів - проекції сторін елемента - постійні величини.

Визначальні співвідношення. Визначальне співвідношення для анізотропного ідеальнопружного матеріалу записується у вигляді добре відомого закону Гука. У матричному вигляді воно було одержане раніше. Тут ми його ще раз приведемо для повноти викладення

(3.85)

де - пружні модулі матеріалу.

Якщо вважати, що пружні властивості матеріалу не змінюються на елементі, то, оскільки деформації постійні, напруження також постійні на елементі.

Елементні матриці жорсткості. Для задачі про плоскодеформований стан тіла загальний вигляд матриць жорсткості був одержаний у попередніх розділах. Нагадаємо його:

. (3.86)

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

. (3.87)

Якщо елемент має постійну товщину (h=const), то h також виноситься за знак інтеграла за площею елемента, який, таким чином, легко обчислюється і буде дорівнювати площі елемента А. У результаті одержимо вираз елементної матриці жорсткості

(3.88)

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

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

. (3.89)

Згідно із загальною формулою обчислення вектора об'ємних вузлових сил

. (3.90)

Нехай товщина елемента стала і густина об'ємних сил також стала. Тоді h і b можуть бути винесені за знак інтеграла і необхідно буде обчислити тільки інтеграли від за площею елемента

. (3.91)

Тут можна скористатися відомою математичною формулою обчислення інтеграла за площею трикутника

. (3.92)

Таким чином, розв'язавши інтеграли від функцій інтерполяції і перемноживши одержану матрицю на вектор b, одержимо остаточний вираз елементного вектора вузлових сил, який статично еквівалентний заданим постійним об'ємним силам

. (3.93)

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

3.9 Ізопараметричний підхід у МСЕ

Вступ. Техніка побудови лінійного трикутного скінченного елемента, викладена в попередніх розділах, формально може бути поширена на чотиристоронні або тристоронні елементи високого порядку. Проте це негайно призведе до ряду труднощів:

1) побудова функцій форми (функцій інтерполяції), які задовольняють всі необхідні умови, для елементів з криволінійними границями є дуже складною задачею;

2) інтеграли, які виникають у виразі елементних матриць жорсткості і вузлових векторів, вже не можуть бути обчислені в замкнутій формі.

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

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

Початкові поняття. Розглянутий у попередніх розділах лінійний трикутний скінченний елемент є простим ізопараметричним елементом, хоча його рівняння і були одержані іншим шляхом. Запишемо основні співвідношення цього елемента:

, (3.94а)

(3.94б)

. (3.94в)

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

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

Повернемося до розгляду ізопараметричного елемента. Перепишемо співвідношення (3.94а-в) в універсальному вигляді, щоб підкреслити їх внутрішню єдність:

(3.95)

Формула (3.95) підкреслює той факт, що геометрія і переміщення задані одним і тим самим параметричним поданням за допомогою функцій форми . У цьому і полягає ключова ідея ізопараметричного підходу в методі скінченних елементів.

Ця методика легко розповсюджується на випадок довільного за геометрією елемента, включаючи і тривимірні елементи. При цьому термін "трикутні координати" замінюється на загальніше поняття "природні локальні координати", які також можуть бути інтерпретовані як параметри, які визначають відповідні функції форми елемента.

Загальне ізопараметричне формулювання пружного двовимірного елемента. Розглянемо загальне ізопараметричне формулювання пружного двовимірного елемента, який складається з n вузлів. Для опису елемента потрібні два типи співвідношень інтерполяції: одне - для задання геометрії елемента, а друге - для апроксимації вектора переміщень на елементі. Обидва співвідношення будуються за допомогою однакових функцій форми і відповідних значень функцій (глобальних координат і компонент переміщень) у вузлах елемента. Таким чином, одержуємо:

1) геометричні співвідношення (інтерполяцію координат):

(3.96)

2) функціональні співвідношення (інтерполяцію переміщень):

. (3.97)

Рівняння (3.96) і (3.97) можуть бути подані у вигляді єдиної формули в матричному вигляді:

(3.98)

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

. (3.99)

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

Лінійний ізопараметричний трикутний елемент. Розглянемо типові ізопараметричні скінченні елементи, які широко використовуються в скінченно-елементних програмах для розв'язання двовимірних задач механіки. Лінійний трикутний елемент може бути поданий тепер як простий ізопараметричний структурний двовимірний елемент, який описується матричним рівнянням

(3.100)

Функції форми виражаються елементарними співвідношеннями через трикутні координати:

. (3.101)

Його вигляд показаний на рис. 3.16.

Рисунок 3.16 - Тривузловий лінійний ізопараметричний трикутний елемент

Квадратичний ізопараметричний трикутний елемент На рисунку 3.17 зображений шестивузловий трикутний елемент, що також є типовим представником сім'ї ізопараметричних скінченних елементів, які широко застосовуються в розрахунках. Він описується таким матричним рівнянням:

. (3.102)

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

(3.103)

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

Рисунок 3.17 - Шестивузловий трикутний скінченний елемент: (а) - суперпараметричний варіант елемента, що має прямолінійні сторони і середні вузли 4, 5, 6, які розміщенні посередині відповідних сторін елемента; (б) - ізопараметричний варіант того самого елемента

3.10 Алгоритм МСЕ для тривимірної задачі теорії пружності

Вступ. Динамічна задача лінійної теорії в'язкопружності (рис. 3.18) описується добре відомою системою тензорних рівнянь в заданому об'ємі суцільного середовища V з граничними умовами кінематичного або силового типу на її поверхні

S = S₁ + S₂:

,

,

, (3.104)

де - тензори напружень і деформацій;

e - девіатор тензора деформацій

;

E - одиничний тензор;

u=u(r,t) - вектор переміщень, який залежить від часу t і радіуса вектора точок суцільного середовища r;

f_v=f_v(r,t) - вектор заданих об'ємних сил, який залежить від часу t і радіуса вектора точок суцільного середовища r;