Метод скінченних елементів у механіці конструкцій

Походження методу кінцевих елементів простежується до матричного аналізу структур[1][2] , де було введено поняття матриці переміщення та матриці жорсткості. Концепції скінченних елементів були розвинені на основі інженерних методів в 1950-х роках. Метод скінченних елементів дуже розвинувася у 1960-х і 1970-х роках завдяки Джону Арджірісі і його колегам(університеті Штутгарта), Раю В. Клафу ( університет Каліфорнії, Берклі), О. Зінкевичу та іншим. Оригінальні твори, такі, як Арджіріс [3] і Клаф [4] і стали основою методів сучасного скінченноелементного структурного аналізу. Такі книги, як Зінкеивич [5] та більш пізні книги, такі як Янг [6] дають докладно розповідають про історію розвитку методу скінченних елементів. Реалізація методу у програмних продуктах описується в класичному тексті Сміт, Гріффітс і Маргеттс.[7]

• Прямі або зігнуті одновимірні елементи з фізичними властивостями, такими як осьова, крутильна жорсткість та жорсткість при згині. Такий тип елементів призначений для моделювання кабелів, розкосів, балок, сіток і каркасів. Прямі елементи, як правило, мають два вузли, по одному на кожному кінці, в той час як вигнуті елементи потребують щонайменше три вузли, включаючи вузли на краях. 
• Двовимірні елементи, використовуються лише для сил на площині (напруга на площині, площина деформації) для мембран, і пластин. Вони можуть бути різної форми, плоскими або криволінійними трикутниками і чотирикутниками. Вузли зазвичай розташовуються в краях елемента, і при необхідності для більшої точності, додаткові вузли можуть розташовуватися вздовж елемента краями або навіть всередині елемента. Елементи розташовуються в серединній поверхні.
• Елементи форми тора для осесиметричних задач для мембран, товстих пластин, оболонок і  твердих тіл. У поперечному перерізі елементи схожі на описані раніше типи: одновимірний для тонких пластин і оболонок, а також двовимірних твердих тіл, товстих плит і оболонок.
• Тривимірні елементи для моделювання 3-D твердих тіл, таких як компоненти машин, дамб, насипів і ґрунтових масивів. Звичною формою елемента є tetrahedrals і hexahedrals. Вузли розташовані у вершинах або на поверхнях чи всередині елемента.

З точки зору програмування важливо змоделювати систему так, що:

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

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

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

Принцип віртуальних переміщень закритої системи виражає тотожність зовнішньої і внутрішньої віртуальної роботи:

${\displaystyle {\mbox{External virtual work}}=\int _{V}\delta {\boldsymbol {\epsilon }}^{T}{\boldsymbol {\sigma }}\,dV\qquad \mathrm {(1)} }$

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

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

${\displaystyle \mathbf {R} =\mathbf {Kr} +\mathbf {R} ^{o}\qquad \qquad \qquad \mathrm {(2)} }$

де

${\displaystyle \mathbf {R} }$ = вектор вузлових сил, що представляють зовнішні сили, прикладені до вузлів системи.

${\displaystyle \mathbf {K} }$ = матриця жорсткості, яка є колективною дією окремих елементів матриці жорсткості:${\displaystyle \mathbf {k} ^{e}}$.

${\displaystyle \mathbf {r} }$ = вектор вузлових переміщень системи.

${\displaystyle \mathbf {R} ^{o}}$ = вектор еквівалентних вузлових сил, які представляють усі зовнішні ефекти, крім вузлових сил, які вже включені. Ці зовнішні ефекти можуть включати розподілені або концентровані поверхневі сили, термічні впливи, початкові напруження і деформації.

Після врахування обмежень, вузлові переміщення знаходяться шляхом розв'язання системи лінійних рівнянь (2), тобто:

${\displaystyle \mathbf {r} =\mathbf {K} ^{-1}(\mathbf {R} -\mathbf {R} ^{o})\qquad \qquad \qquad \mathrm {(3)} }$

Потім, деформації і напруження на окремих елементах можуть бути знайдені таким чином:

${\displaystyle \mathbf {\epsilon } =\mathbf {Bq} \qquad \qquad \qquad \qquad \mathrm {(4)} }$

${\displaystyle \mathbf {\sigma } =\mathbf {E} (\mathbf {\epsilon } -\mathbf {\epsilon } ^{o})+\mathbf {\sigma } ^{o}=\mathbf {E} (\mathbf {Bq} -\mathbf {\epsilon } ^{o})+\mathbf {\sigma } ^{o}\qquad \qquad \qquad \mathrm {(5)} }$

де

${\displaystyle \mathbf {q} }$ = вектор вузлових переміщень--підмножина системи переміщення вектора r , що стосується розглянутих елементів.

${\displaystyle \mathbf {B} }$ = матриця деформацій і переміщення, яка перетворює вузлові переміщення  q у напруження в будь-якій точці елемента.

${\displaystyle \mathbf {E} }$ = матриця пружності.

${\displaystyle \mathbf {\epsilon } ^{o}}$ = вектор початкових деформацій в елементах.

${\displaystyle \mathbf {\sigma } ^{o}}$ = вектор початкових напружень в елементах.

Застосовуючи рівняння віртуальної роботи (1) до системи, отримуємо матриці елементів ${\displaystyle \mathbf {B} }$, ${\displaystyle \mathbf {k} ^{e}}$ а також схему побудови системи матриць ${\displaystyle \mathbf {R} ^{o}}$ і ${\displaystyle \mathbf {K} }$. Інші матриці, такі як  ${\displaystyle \mathbf {\epsilon } ^{o}}$, ${\displaystyle \mathbf {\sigma } ^{o}}$, ${\displaystyle \mathbf {R} }$ and ${\displaystyle \mathbf {E} }$ є відомими і можуть бути безпосередньо створені з вхідних даних.

Нехай ${\displaystyle \mathbf {q} }$ - вектор вузлових переміщень стандартного елемента. Переміщення в будь-якій іншій точці елемента може бути знайдений за рахунок використання інтерполяції функцій:

${\displaystyle \mathbf {u} =\mathbf {N} \mathbf {q} \qquad \qquad \qquad \mathrm {(6)} }$

де

${\displaystyle \mathbf {u} }$ = вектора переміщень в будь-якій точці {х,у,z} елемента.

${\displaystyle \mathbf {N} }$ = матриця базисних функцій, яка виступає інтерполяційною функцією.

Рівняння (6) утворює нові величини:

• Віртуальні переміщення, які є функцією віртуальних вузлових переміщень: ${\displaystyle \delta \mathbf {u} =\mathbf {N} \delta \mathbf {q} \qquad \qquad \qquad \mathrm {(6b)} }$
• Деформації в елементах, що виникають в результаті переміщення вузлів елемента: ${\displaystyle \mathbf {\epsilon } =\mathbf {Du} =\mathbf {DNq} \qquad \qquad \qquad \qquad \mathrm {(7)} }$

де ${\displaystyle \mathbf {D} }$ = матриця диференціальних операторів перетворення переміщень в деформації з допомогою лінійної теорії пружності теорії. Рівняння (7) показує, що матриця B у (4) обчислюється таким чином:

${\displaystyle \mathbf {B} =\mathbf {DN} \qquad \qquad \qquad \qquad \mathrm {(8)} }$

• Віртуальні деформації утворені, віртуальними переміщення вузлових елементів: ${\displaystyle \delta {\boldsymbol {\epsilon }}=\mathbf {B} \delta \mathbf {q} \qquad \qquad \qquad \qquad \mathrm {(9)} }$

Для стандартного елемента об'єму ${\displaystyle V^{e}}$, внутрішня віртуальна робота через віртуальних переміщення знаходиться шляхом підстановки (5) і (9) в (1):

${\displaystyle {\mbox{Internal virtual work}}=\int _{V^{e}}\delta {\boldsymbol {\epsilon }}^{T}{\boldsymbol {\sigma }}\,dV^{e}=\delta \ \mathbf {q} ^{T}\int _{V^{e}}\mathbf {B} ^{T}{\big \{}\mathbf {E} (\mathbf {Bq} -\mathbf {\epsilon } ^{o})+\mathbf {\sigma } ^{o}{\big \}}\,dV^{e}\qquad \mathrm {(10)} }$

Сумуючи внутрішню віртуальну роботу (14) для всіх елементів дає праву частину (1):

${\displaystyle {\mbox{System internal virtual work}}=\sum _{e}\delta \ \mathbf {r} ^{T}{\big (}\mathbf {k} ^{e}\mathbf {r} +\mathbf {Q} ^{oe}{\big )}=\delta \ \mathbf {r} ^{T}{\big (}\sum _{e}\mathbf {k} ^{e}{\big )}\mathbf {r} +\delta \ \mathbf {r} ^{T}\sum _{e}\mathbf {Q} ^{oe}\qquad \mathrm {(15)} }$

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

• Робота вузлових сил R: ${\displaystyle \delta \ \mathbf {r} ^{T}\mathbf {R} \qquad \mathrm {(16)} }$
• Робота зовнішніх сил ${\displaystyle \mathbf {T} ^{e}}$ на частині ${\displaystyle \mathbf {S} ^{e}}$ ребер або поверхонь елементів і сили ${\displaystyle \mathbf {f} ^{e}}$

${\displaystyle \sum _{e}\int _{S^{e}}\delta \ \mathbf {u} ^{T}\mathbf {T} ^{e}\,dS^{e}+\sum _{e}\int _{V^{e}}\delta \ \mathbf {u} ^{T}\mathbf {f} ^{e}\,dV^{e}}$

Підстановляючи у (6b) отримуємо:

${\displaystyle \delta \ \mathbf {q} ^{T}\sum _{e}\int _{S^{e}}\mathbf {N} ^{T}\mathbf {T} ^{e}\,dS^{e}+\delta \ \mathbf {q} ^{T}\sum _{e}\int _{V^{e}}\mathbf {N} ^{T}\mathbf {f} ^{e}\,dV^{e}}$

або ${\displaystyle -\delta \ \mathbf {q} ^{T}\sum _{e}(\mathbf {Q} ^{te}+\mathbf {Q} ^{fe})\qquad \mathrm {(17a)} }$

де матриці для елементів визначені як:

${\displaystyle \mathbf {Q} ^{te}=-\int _{S^{e}}\mathbf {N} ^{T}\mathbf {T} ^{e}\,dS^{e}\qquad \mathrm {(18a)} }$

${\displaystyle \mathbf {Q} ^{fe}=-\int _{V^{e}}\mathbf {N} ^{T}\mathbf {f} ^{e}\,dV^{e}\qquad \mathrm {(18b)} }$

Знову ж таки, чисельне інтегрування часто застосовується для їх обчислення. Аналогічно, заміняючи q у (17a) на r , після перетворення:

${\displaystyle -\delta \ \mathbf {r} ^{T}\sum _{e}(\mathbf {Q} ^{te}+\mathbf {Q} ^{fe})\qquad \mathrm {(17b)} }$

Додавання (16), (17b) і прирівнювання суми до (15) дає: ${\displaystyle \delta \ \mathbf {r} ^{T}\mathbf {R} -\delta \ \mathbf {r} ^{T}\sum _{e}(\mathbf {Q} ^{te}+\mathbf {Q} ^{fe})=\delta \ \mathbf {r} ^{T}{\big (}\sum _{e}\mathbf {k} ^{e}{\big )}\mathbf {r} +\delta \ \mathbf {r} ^{T}\sum _{e}\mathbf {Q} ^{oe}}$

Так як віртуальні переміщення ${\displaystyle }$ є довільними, рівність зводиться до:

${\displaystyle \mathbf {R} ={\big (}\sum _{e}\mathbf {k} ^{e}{\big )}\mathbf {r} +\sum _{e}{\big (}\mathbf {Q} ^{oe}+\mathbf {Q} ^{te}+\mathbf {Q} ^{fe}{\big )}}$

Порівняння з (2) показує, що:

• Матриця жорсткості системи отримується  шляхом підсумовування матриці жорсткості для кожного елемента:

${\displaystyle \mathbf {K} =\sum _{e}\mathbf {k} ^{e}}$

• Вектор еквівалентних сил для вузлів отримується шляхом підсумовування елементів векторів навантаження:

${\displaystyle \mathbf {R} ^{o}=\sum _{e}{\big (}\mathbf {Q} ^{oe}+\mathbf {Q} ^{te}+\mathbf {Q} ^{fe}{\big )}}$

На практиці, матриці елементів ніколи не розписують і не розширюють. Навпаки, матриця жорсткості для системи  ${\displaystyle \mathbf {K} }$ отримується додаванням коефіцієнтів ${\displaystyle {k}_{ij}^{e}}$ to ${\displaystyle {K}_{kl}}$ , де індекси ij, kl означають, що переміщення даного елемента  ${\displaystyle {q}_{i}^{e},{q}_{j}^{e}}$ збігаються, відповідно, з переміщенням вузлів системи ${\displaystyle {r}_{k},{r}_{l}}$. Подібно, ${\displaystyle \mathbf {R} ^{o}}$ ${\displaystyle {Q}_{i}^{e}}$ до ${\displaystyle {R}_{k}^{o}}$ , де ${\displaystyle {q}_{i}^{e}}$ відповідає ${\displaystyle {r}_{k}}$. Пряме додавання ${\displaystyle {k}_{ij}^{e}}$ до ${\displaystyle {K}_{kl}}$ породжує процедуру  Direct Stiffness Method.

• Львівський національний університет ім. І. Франка