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

1. Метод прямої жорсткості

Будівельний МСЕ з'явився за десятиліття до загального континуального МСЕ — він виріс із матричного аналізу каркасних конструкцій, розробленого для авіаційного проєктування у 1950-х роках. Основна ідея проста: моделювати конструкцію як мережу одновимірних елементів (стрижнів, балок, колон), з'єднаних у дискретних вузлах. Кожен елемент має відоме співвідношення між силами на його кінцях і переміщеннями на цих кінцях — його матрицю жорсткості елемента. Об'єднавши матриці жорсткості всіх елементів в одну велику глобальну матрицю жорсткості, прикладаємо навантаження й опори та розв'язуємо лінійну систему для невідомих переміщень.

Це та сама система K·u = F, яку розв'язує континуальний МСЕ, але тут k_e отримано з точного аналітичного розв'язку керівного ОДР для стрижня чи балки, а не шляхом чисельного інтегрування функцій форми по довільній формі елемента. Саме тому ручний розрахунок невеликих ферм був практичним ще до появи комп'ютерів — і саме тому метод прямої жорсткості й досі є найшвидшим та найпрозорішим способом навчання основам МСЕ.

2. Елементи ферми (стрижні)

Стрижень ферми сприймає лише осьову силу — розтяг або стиск уздовж своєї довжини — і шарнірно закріплений з обох кінців, тому не передає згинального моменту. У власній локальній системі координат (x′ уздовж стрижня) двовузловий стрижень з осьовою жорсткістю EA/L має добре відому одновимірну матрицю жорсткості:

У 2D або 3D кожен вузол має 2 або 3 поступальні ступені свободи, тому локальну осьову жорсткість потрібно повернути у глобальну систему координат x-y(-z), використовуючи напрямні косинуси стрижня. Ферменна модель мосту чи даху може мати кілька сотень стрижнів і вузлів — достатньо мало, щоб розраховувати вручну в докомп'ютерну епоху, саме тому метод спочатку розробили для ферм.

3. Елементи балки Ейлера-Бернуллі

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

Ермітові кубічні функції форми, використані тут, інтерполюють одночасно переміщення та нахил, гарантуючи неперервність обох на межах елементів (C¹-неперервність) — це необхідно, оскільки керівне рівняння EI·v'''' = q є рівнянням четвертого порядку. Елементи балки Тимошенка розширюють цю модель, додаючи деформацію поперечного зсуву — це важливо для коротких, глибоких балок, де гнучкість на зсув суттєва (відношення прольоту до висоти нижче ~10).

4. Елементи рами та перетворення координат

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

У 3D елементи рами також несуть кручення (GJ/L, де G — модуль зсуву, а J — момент інерції при крученні) та згин відносно двох ортогональних осей, що дає повну матрицю жорсткості елемента 12×12 (6 СС на вузол: 3 поступальні, 3 обертальні). Це базовий елемент, на якому працює майже кожен пакет для розрахунку конструкцій — SAP2000, ETABS, STAAD.Pro та Robot Structural Analysis будують свою глобальну матрицю жорсткості саме з такого 12×12 елемента рами.

5. Збірка та розв'язання глобальної системи

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

K розріджена та стрічкова — вузол з'єднується лише з кількома сусідніми елементами, тому більшість записів дорівнює нулю. Перенумерація СС для мінімізації ширини стрічки (Cuthill-McKee, вкладена диссекція) значно прискорює прямі розв'язувачі. Коли D_f відоме, локальні переміщення кожного елемента відновлюються, а внутрішні епюри осьової сили, поперечної сили та згинального моменту випливають безпосередньо з k_local · d_e.

6. Статична конденсація та підструктуризація

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

Це той самий метод редукції Гайана (Guyan), що використовується для стиснення величезних скінченно-елементних моделей у динамічному аналізі, і саме так інтегроване з САПР програмне забезпечення для рам збирає попередньо обчислені "суперелементи" для типових компонентів (збірна балка, опора мосту) без повторного виведення їхньої внутрішньої поведінки щоразу.

7. Модальний аналіз і власні частоти

Конструкції не просто несуть статичні навантаження — вони вібрують під дією вітру, сейсмічного збудження, кроків людей чи роботи обладнання. Додавання узгодженої матриці мас M перетворює статичну задачу на узагальнену задачу на власні значення:

Розв'язання для найменших власних значень дає основні власні частоти та форми коливань, що домінують у динамічній відповіді конструкції. Резонанс виникає, коли частота зовнішнього збудження наближається до ω — саме тому висотні будівлі перевіряють на частоти вихрового зриву від вітру, а пішохідні мости (сумнозвісний Millennium Bridge у Лондоні у 2000 році) перевіряють на частоти кроків пішоходів близько 2 Гц.

8. Лінійний аналіз стійкості (втрата стійкості)

Стрункі стиснуті елементи можуть зазнавати втрати стійкості — раптової бокової нестійкості — при навантаженнях, значно нижчих за межу текучості матеріалу. Лінійний (власне-значеневий) аналіз стійкості вводить геометричну матрицю жорсткості K_g, яка враховує, як осьове навантаження знижує ефективну згинальну жорсткість:

Аналіз стійкості методом МСЕ узагальнює формулу Ейлера для колони на довільні рами, де нестійкість може проявлятися як зсув цілого поверху, а не згин окремого елемента. Найменше власне значення λ₁ дає коефіцієнт критичного навантаження; будівельні норми зазвичай вимагають, щоб λ₁ мало суттєвий запас понад 1.0, а нелінійний (великодеформаційний) аналіз використовують, коли важлива поведінка після втрати стійкості або чутливість до недосконалостей.

9. Реалізація 2D розв'язувача рам

Мінімальному 2D розв'язувачу жорсткості рам потрібні лише три складові: генерація жорсткості елемента, збірка та обробка граничних умов. Ось ядро процедури обчислення жорсткості елемента на JavaScript, що відповідає матриці 6×6, виведеній у Розділі 4:

function frameElementStiffness(E, A, I, L, angleRad) {
  // локальна жорсткість 6x6: [u1,v1,th1,u2,v2,th2]
  const k = zeros(6, 6);
  const EA_L  = E * A / L;
  const EI_L3 = E * I / (L ** 3);
  const EI_L2 = E * I / (L ** 2);
  const EI_L  = E * I / L;

  k[0][0] =  EA_L;  k[3][3] =  EA_L;  k[0][3] = -EA_L;  k[3][0] = -EA_L;

  k[1][1] =  12 * EI_L3;  k[4][4] =  12 * EI_L3;
  k[1][4] = -12 * EI_L3;  k[4][1] = -12 * EI_L3;

  k[1][2] = k[2][1] =  6 * EI_L2;
  k[1][5] = k[5][1] =  6 * EI_L2;
  k[2][4] = k[4][2] = -6 * EI_L2;
  k[4][5] = k[5][4] = -6 * EI_L2;

  k[2][2] = 4 * EI_L;  k[5][5] = 4 * EI_L;
  k[2][5] = k[5][2] = 2 * EI_L;

  const c = Math.cos(angleRad), s = Math.sin(angleRad);
  const T = buildTransform(c, s); // блок-діагональна матриця повороту 6x6
  return multiply(transpose(T), multiply(k, T)); // Tᵀ k T
}

function assembleGlobal(elements, nNodes) {
  const ndof = nNodes * 3;
  const K = zeros(ndof, ndof);
  for (const el of elements) {
    const ke = frameElementStiffness(el.E, el.A, el.I, el.L, el.angle);
    const map = [el.n1*3, el.n1*3+1, el.n1*3+2,
                el.n2*3, el.n2*3+1, el.n2*3+2];
    for (let i = 0; i < 6; i++)
      for (let j = 0; j < 6; j++)
        K[map[i]][map[j]] += ke[i][j];
  }
  return K; // далі розбити на вільні/опорні СС і розв'язати K_ff Df = Ff
}

Далі веб-візуалізатору потрібні ще три кроки: (1) розв'язати K_ff·D_f = F_f прямим або розрідженим лінійним розв'язувачем (метод Гаусса цілком підходить для кількох тисяч СС), (2) відновити кінцеві зусилля елементів з k_local·d_e для кожного елемента, і (3) відобразити деформовану форму та епюри згинальних моментів, масштабуючи та інтерполюючи ермітову кубічну функцію всередині кожного елемента балки — саме такий інтерактивний робочий процес використовується у симуляціях конструкцій на цьому сайті.