Метод скінченних елементів (МСЕ) простими словами

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

1. Головна ідея

Більшість фізичних задач зводиться до диференціальних рівнянь у частинних похідних (ДРЧП): напруження в балці, тепло крізь стіну, обтікання крила повітрям. Ці рівняння рідко мають аналітичні розв’язки для складних геометрій. МСЕ розв’язує їх чисельно, виконуючи:

Дискретизація

Розбити область на маленькі елементи (трикутники, тетраедри)

Апроксимація

Припустити, що розв’язок змінюється поліноміально в межах кожного елемента

Збирання

Об’єднати всі елементи в глобальну систему рівнянь

Розв’язання

За допомогою лінійної алгебри знайти невідомі значення в кожному вузлі

Більше елементів означає кращу точність, але більше обчислень. Сучасна автомобільна краш-симуляція використовує 5–10 мільйонів елементів і триває години на кластері GPU.

2. Сильна форма → слабка форма

Розгляньмо 1D стрижень під осьовим навантаженням. Сильна (диференціальна) форма рівноваги така:

-d/dx [EA \cdot du/dx] = f(x) для x \in [0, L] E = модуль Юнга, A = площа поперечного перерізу, u = переміщення, f = розподілене навантаження

Сильна форма вимагає, щоб u було двічі диференційовним усюди — занадто суворо для чисельних методів. Множимо на тестову функцію v(x) та інтегруємо частинами, щоб отримати слабку форму:

\int₀ᴸ EA \cdot (du/dx) \cdot (dv/dx) dx = \int₀ᴸ f\cdotv dx + [EA\cdot(du/dx)\cdotv]₀ᴸ Це знижує вимогу до гладкості: u та v потребують лише перших похідних (неперервність C⁰), а не других.

Слабка форма — математична основа МСЕ. Вона каже: знайти u таке, щоб рівняння виконувалося для всіх допустимих тестових функцій v. У МСЕ і u, і v апроксимуються одними й тими самими поліноміальними функціями форми — це метод Гальоркіна.

3. Побудова сітки області

Область (балка, лопатка турбіни, череп) поділяється на елементи, що не перекриваються:

1D: відрізки прямих (стрижні, балки)
2D: трикутники або чотирикутники (пластини, оболонки)
3D: тетраедри, гексаедри, призми або піраміди

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

Збіжність сітки: у міру зменшення елементів (h → 0) розв’язок МСЕ збігається до точного розв’язку ДРЧП. Швидкість збіжності залежить від типу елемента: лінійні елементи збігаються як O(h²) за переміщенням, квадратичні — як O(h³). Завжди проводьте дослідження збіжності сітки, перш ніж довіряти результатам.

4. Функції форми та інтерполяція

У межах кожного елемента невідоме поле u апроксимується зваженою сумою функцій форми N_i:

u(x) = Σᵢ Nᵢ(x) \cdot uᵢ де uᵢ — вузлові значення (невідомі) а Nᵢ(x) — поліноміальні функції форми: - Лінійна (2 вузли): N₁ = (1-ξ)/2, N₂ = (1+ξ)/2 (ξ \in [-1,1]) - Квадратична (3 вузли): додає серединний вузол - Ізопараметрична: ті самі функції відображають геометрію та поле

Функції форми мають властивість розкладу одиниці: Σ Nᵢ = 1 усюди. Кожна Nᵢ дорівнює 1 у власному вузлі та 0 в усіх інших вузлах. Це забезпечує неперервність на межах елементів (C⁰ для елементів Лагранжа).

5. Матриця жорсткості

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

kᵉ \cdot uᵉ = fᵉ Матриця жорсткості елемента: kᵉᵢⱼ = \int EA \cdot (dNᵢ/dx) \cdot (dNⱼ/dx) dx Вектор сил елемента: fᵉᵢ = \int f \cdot Nᵢ dx + граничні доданки

Для 2-вузлового 1D лінійного елемента довжиною Lᵉ:

kᵉ = (EA/Lᵉ) \cdot [ 1 -1 ] [-1 1 ] Це матриця «жорсткості пружини» — базовий будівельний блок.

Для 2D та 3D елементів інтеграли обчислюються чисельно за допомогою квадратури Гаусса — підсумовуванням підінтегральної функції в певних точках вибірки з певними вагами. Чотирикутний елемент зазвичай використовує 2×2 точки Гаусса; гексаедр використовує 2×2×2 = 8 точок Гаусса.

6. Збирання та розв’язання

Глобальна матриця жорсткості K збирається відображенням локальної матриці кожного елемента в глобальну нумерацію ступенів свободи за допомогою таблиці зв’язності. У результаті отримуємо велику, розріджену, симетричну додатно визначену систему:

K \cdot U = F K = глобальна матриця жорсткості N\timesN (N = загальна кількість ступенів свободи) U = вектор усіх невідомих вузлових переміщень F = вектор прикладених сил та граничних умов Для мільйона невідомих: K має ~10⁶ рядків, але лише ~50 ненульових на рядок \to зберігання розрідженої матриці, ітеративні розв’язувачі (CG, GMRES)

Після застосування граничних умов (наприклад, u=0 на опорах) система розв’язується. Постобробка обчислює похідні величини: напруження σ = E·ε, деформацію ε = du/dx, еквівалентне напруження за Мізесом для прогнозування руйнування.

7. МСЕ у реальному світі

Структурний аналіз: кожен кузов автомобіля, фюзеляж літака та міст аналізуються методом МСЕ на статичні навантаження, втому та краш-удар. ANSYS, Abaqus, NASTRAN є галузевими стандартами.
Теплопередача: МСЕ розв’язує рівняння теплопровідності, щоб прогнозувати теплові градієнти в електроніці, печах та теплозахисних екранах при поверненні в атмосферу.
Гідродинаміка (CFD): варіанти МСЕ (SUPG, DG) розв’язують рівняння Нав’є–Стокса для складних геометрій. Широко використовуються OpenFOAM (метод скінченних об’ємів) та COMSOL (МСЕ).
Електромагнетизм: МСЕ розв’язує рівняння Максвелла для проєктування антен, оптимізації двигунів та проєктування котушок МРТ.
Біомедицина: аналіз напружень у кістках конкретного пацієнта за КТ-знімками (сітка генерується безпосередньо з воксельних даних), симуляція серцевих клапанів, проєктування зубних імплантів.

МСЕ у браузері: бібліотеки JavaScript, як-от FEA.js та Three.js, можуть візуалізувати сітки та результати МСЕ. WebAssembly робить реальним запуск невеликих розв’язувачів МСЕ повністю на боці клієнта — наша власна симуляція аналізу напружень використовує цей підхід.