🔧 Інженерія · Чисельні методи
📅 Березень 2026⏱ 14 хв🔴 Просунутий
·
Останнє оновлення: 28 травня 2026 р.
Метод скінченних елементів зрозуміло: побудова сітки, розв'язання та результати
Кожен сучасний літак, симуляція автомобільної аварії та конструкторський розрахунок залежать від МСЕ. Він перетворює диференціальні рівняння у частинних похідних, що описують неперервні фізичні області, на системи алгебраїчних рівнянь, які можуть розв'язувати комп'ютери, — поділяючи область на малі елементи й наближаючи невідоме поле в межах кожного з них.
1. Від ДРЧП до алгебраїчної системи
Більшість інженерних задач описуються диференціальними рівняннями у частинних похідних (ДРЧП). Для лінійної пружності:
Визначальне ДРЧП (рівновага, без об'ємних сил):
∇·σ = 0 (дивергенція тензора напружень = 0)
σ = C : ε (визначальне співвідношення: напруження = тензор пружності × деформація)
ε = ½(∇u + (∇u)ᵀ) (співвідношення деформація-переміщення)
Де:
u = поле переміщень (те, що ми шукаємо)
σ = тензор напружень (6 незалежних компонент у 3D)
ε = тензор деформацій
C = тензор пружності 4-го порядку (зводиться до 2 сталих E, ν для ізотропних матеріалів)
Аналітичні розв'язки існують лише для простих геометрій.
МСЕ розв'язує слабку (варіаційну) форму:
∫_Ω δε : σ dV = ∫_Γ δu · t dA + ∫_Ω δu · b dV
δu = віртуальне переміщення (пробна функція)
t = поверхневі навантаження (сили на одиницю площі)
b = об'ємні сили (гравітація тощо)
2. Побудова сітки та її якість
Область поділяється на елементи, що не перекриваються. Типи елементів у 2D: tria3/tria6 (трикутники з 3 або 6 вузлами), quad4/quad8 (чотирикутники). У 3D: tet4/tet10 (тетраедри), hex8/hex20 (гексаедри), wedge6/wedge15 (призми).
Ключові показники якості:
- Співвідношення сторін: Відношення найдовшого ребра до найкоротшого. Велике співвідношення сторін (дуже видовжені елементи) погіршує точність. Ціль < 5 для більшості елементів, < 10 лише у примежових шарах.
- Відношення якобіана: Відношення максимального до мінімального визначника якобіана в межах елемента. Має бути > 0 усюди (без вивернутих елементів). Зниження до менш ніж 0.6 спричиняє низьку точність.
- Кутова якість: Для трикутників/тетраедрів дуже малі чи дуже великі кути знижують точність. Ідеальний трикутник: рівносторонній. Тріангуляція Делоне максимізує мінімальний кут.
- Перехід розмірів елементів: Поступові переходи розміру (коефіцієнт росту < 1.5) між дрібними областями (напр., концентраторами напружень) і грубими областями.
Сітка примежового шару (CFD): У динаміці рідин швидкість змінюється від нуля (стінка) до швидкості незбуреного потоку в межах тонкого примежового шару. Це вимагає призматичних шарів дуже тонких елементів з великим співвідношенням сторін біля стінок — зазвичай 20+ шарів з y+ ≈ 1 для прямого розрахунку біля стінки в моделюванні турбулентності. Автоматична генерація таких сіток зі складної CAD-геометрії — одне з найскладніших завдань побудови сітки.
3. Функції форми
У межах кожного елемента невідоме поле u наближається як зважена сума функцій форми N_i(x), де вагами є вузлові значення u_i:
u(x) ≈ Σᵢ Nᵢ(x) · uᵢ (у межах одного елемента)
Функції форми для лінійного трикутника (CST — стала деформація):
L₁, L₂, L₃ = площинні (барицентричні) координати
N₁ = L₁, N₂ = L₂, N₃ = L₃
Властивості:
Nᵢ(xⱼ) = δᵢⱼ (1 у вузлі i, 0 в усіх інших вузлах)
Σᵢ Nᵢ = 1 (розбиття одиниці)
Функції форми серендипного чотирикутника (Q8) у кутовому вузлі:
N₁(ξ,η) = ¼(1+ξ₀)(1+η₀)(ξ₀+η₀−1)
де ξ₀ = ξ·ξ₁, η₀ = η·η₁ (натуральні координати від −1 до +1)
Функції форми вищого порядку (p-уточнення):
Збільшення степеня полінома підвищує точність без зменшення елементів.
p-МСЕ може досягати експоненційної збіжності для гладких задач.
4. Жорсткість елемента та складання
Для кожного елемента матриця жорсткості елемента K_e обчислюється підстановкою функцій форми у слабку форму:
Матриця жорсткості елемента:
K_e = ∫_Ωe Bᵀ C B dV
B = матриця деформація-переміщення = ∂N (похідні функцій форми)
C = матриця пружності (сталі E, ν)
Чисельне інтегрування (квадратура Гаусса):
∫_Ωe f dV ≈ Σₖ wₖ · f(xₖ) · det(J(xₖ))
xₖ = точки Гаусса, wₖ = ваги, J = якобіан відображення
Повна система МСЕ:
K · u = F
K = глобальна матриця жорсткості (складена з Kₑ)
u = вектор вузлових переміщень (невідомі)
F = вектор сил (прикладені навантаження)
Для моделі з N вузлами та 3 ступенями вільності на вузол:
Розмір матриці: 3N × 3N
Типова сітка: 1,000–10,000,000 вузлів
K розріджена (ненульові лише зв'язані елементи)
Для великих задач потрібні розв'язувачі розріджених систем (PARDISO, MUMPS)
5. Граничні умови та розв'язання
Граничні умови задають відомі величини на межі області. Для міцнісного аналізу:
- Умова Діріхле (головна): Задані переміщення (напр., нерухома опора: u = 0). Застосовується модифікацією рядків/стовпців матриці жорсткості.
- Умова Неймана (природна): Задані сили/навантаження. З'являються природно у векторі сил F.
- Умови контакту: Обмеження-нерівності (поверхні не можуть взаємопроникати). Потребують ітераційних методів розв'язання (штрафний метод, множники Лагранжа або mortar-контакт).
Для лінійних задач K·u = F розв'язується один раз прямими (LU-розклад, Холецький) або ітераційними (метод спряжених градієнтів, GMRES) методами. Для нелінійних задач (великі деформації, пластичність, контакт) ітерація Ньютона–Рафсона розв'язує повторно, доки нев'язка < допуску.
6. Збіжність та оцінювання похибки
МСЕ дає наближені розв'язки. Точність зростає з:
- h-уточнення: Зменшення розміру елемента h. Для лінійних елементів похибка переміщень ∝ h², похибка напружень ∝ h.
- p-уточнення: Збільшення порядку полінома p. Експоненційна збіжність для гладких задач.
- hp-уточнення: Поєднання обох. Найефективніше для задач із локальними сингулярностями (вершини тріщин, вхідні кути).
Апостеріорна оцінка похибки (Зенкевич–Чжу):
η_e = || σ_h* − σ_h ||_Ωe / || σ_h ||_Ω
σ_h* = згладжене напруження (відновлене за SPR або PPR)
σ_h = «сире» напруження МСЕ (розривне на межах елементів)
η = глобальна оцінка похибки (сума за елементами)
Ціль: η < 5-10% для інженерної точності
Використовуйте оцінку похибки для адаптивного уточнення сітки:
Уточнюйте елементи, де локальна похибка > порога
Згрублюйте елементи, де локальна похибка << порога
7. Мультифізичні застосування
- Симуляція аварії (нелінійний нестаціонарний МСЕ): Автомобільна аварія триває ~100 ms. Явне інтегрування за часом із кроками ~1 μs. 10–50 мільйонів елементів. Враховує великі деформації, пластичність, руйнування, контакт та наповнення подушки безпеки. Інструменти: Abaqus/Explicit, LS-DYNA, PAM-CRASH.
- Втома конструкцій: МСЕ дає поля напружень для аналізу за кривою S-N та механіки руйнування. Підрахунок циклів методом «дощового потоку» + правило Майнера оцінюють ресурс деталі за навантаження змінної амплітуди.
- Термомеханічне зв'язування: МСЕ розв'язує теплопередачу (теплопровідність, конвекцію), щоб отримати поле температур, а потім використовує теплові деформації у міцнісному аналізі. Критично для проєктування лопаток турбін (Δ T = 800°C → теплові напруження домінують над механічним навантаженням).
- Електромагнітний МСЕ: Розв'язує рівняння Максвелла у частотній області. Застосовується для проєктування антен, характеристик двигунів, втрат у осерді трансформатора та проєктування котушок МРТ. Інструменти: COMSOL, Ansys HFSS/Maxwell.
- Біомедицина: Індивідуальний ризик перелому кістки за КТ-знімками, механічне проєктування стентів, контактні напруження в колінному імпланті, взаємодія рідина-структура в кохлеарному імпланті.
- Геомеханіка: Ущільнення колектора, стійкість розломів за зміни порового тиску, проєктування обробки тунелів, безпека гребель. Зв'язаний гідромеханічний МСЕ моделює течію рідини у пористих деформівних середовищах.