Матрицы жёсткости и масс

От непрерывного элемента к узловым матрицам

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

u(x) = N(x) q_e

Именно выбор N(x) задаёт закон распределения перемещений внутри элемента. Матрицы возникают после подстановки этой аппроксимации в выражения для энергии деформации и кинетической энергии:

K_e = ∫_V B^TDB dV,    M_e = ∫_V ρN^TN dV,    ε = Bq_e

Для двухузлового стержня типа 112 осевое перемещение интерполируется линейно. При ξ = x/L:

u(x) = (1 − ξ)u′₁ + ξu′₂

Поэтому осевая деформация ε = du/dx постоянна внутри такого элемента. Здесь оператор B содержит первую производную функций формы. Для рамы типа 113 осевое перемещение также линейно, но поперечный прогиб описывается кубическими функциями Эрмита:

v(x) = H₁(x)v′₁ + H₂(x)θ₁ + H₃(x)v′₂ + H₄(x)θ₂

В изгибной части оператор B содержит уже вторые производные функций H. Поэтому кривизна κ = d²v/dx² и, при постоянном EI, изгибающий момент изменяются по длине линейно.

Таким образом, две точки не воспроизводят произвольную деформированную форму: они задают выбранную аппроксимацию. Разбиение конструкции на большее число элементов добавляет узлы и позволяет точнее описать распределение перемещений и деформаций.

Плоский стержень фермы, тип 112

В локальных осях x′, y′ элемент сопротивляется только изменению расстояния между узлами. В реализации сохранены обе поступательные компоненты каждого узла, поэтому локальный вектор имеет четыре компоненты.

k = EA/L

K_local = k · [ 1  0 -1  0
                0  0  0  0
               -1  0  1  0
                0  0  0  0 ]

Согласованная масса распределяет поступательную инерцию между обоими узлами:

M_local = ρAL/6 · [ 2  0  1  0
                    0  2  0  1
                    1  0  2  0
                    0  1  0  2 ]

Коэффициенты 2 и 1 получаются из интеграла ∫N^TN dx. Недиагональные элементы матрицы — не искусственная связь, а следствие того, что движение любой внутренней точки одновременно зависит от перемещений обоих узлов. Поэтому такая масса называется согласованной с выбранными функциями формы.

У свободного стержня ранг K равен 1: лишь одна линейная комбинация перемещений создаёт осевую деформацию. Остальные три независимые формы — два поступательных перемещения и бесконечно малый поворот элемента как единого целого — в линейной постановке не вызывают осевого удлинения. Поэтому упругая энергия для них равна нулю.

Плоская рама Эйлера–Бернулли, тип 113

Рамный элемент объединяет осевую жёсткость и изгиб. Принято допущение Эйлера–Бернулли: поперечный сдвиг не деформирует сечение, поэтому поворот связан с производной прогиба.

Интегрирование производных этих функций формы даёт характерные изгибные коэффициенты 12, 6, 4 и 2 в матрице жёсткости:

a = EA/L,   b = 12EI/L³,   c = 6EI/L²,   d = 4EI/L,   e = 2EI/L

K_local = [ a   0   0  -a   0   0
            0   b   c   0  -b   c
            0   c   d   0  -c   e
           -a   0   0   a   0   0
            0  -b  -c   0   b  -c
            0   c   e   0  -c   d ]

В createStructuralElement.m используется стандартная согласованная матрица масс:

M_local = ρAL/420 ·
[ 140    0       0     70    0        0
    0   156     22L     0   54      -13L
    0   22L     4L²     0   13L     -3L²
   70    0       0    140    0        0
    0    54     13L     0  156      -22L
    0  -13L    -3L²     0  -22L      4L² ]

Та же кубическая интерполяция прогиба используется в кинетической энергии и приводит к коэффициентам 156, 54, 22L, 13L и другим элементам согласованной матрицы масс.

У свободного рамного элемента три движения как твёрдого тела: два переноса и один поворот. Поэтому ранг матрицы жёсткости 6×6 равен 3.

Единицы и смысл коэффициентов

Матрицы с вращательными СС не имеют одной общей физической единицы для всех записей. Например, связь силы с перемещением имеет Н/м, а связь момента с углом — Н·м. Проверять нужно каждую пару «обобщённая сила — обобщённое перемещение».

Пять проверок до сборки

1. Размер: 4×4 для типа 112, 6×6 для типа 113.
2. Конечность: ни один коэффициент не равен NaN или Inf.
3. Симметрия: K = K^T и M = M^T.
4. Ранг: 1 у свободного стержня, 3 у свободной рамы.
5. Масса: при единичном переносе квадратичная форма r^TMr равна ρAL.

Матрица жёсткости свободного элемента обязана быть вырожденной. Это не ошибка: элемент может двигаться без деформации. После сборки и задания опор глобальная матрица должна стать невырожденной; если этого не произошло, в модели осталось незакреплённое движение или внутренний механизм.

Лабораторная 01. Жёсткость и масса наклонного элемента

Цель — связать коэффициенты матриц с геометрией и объяснить, почему осевая жёсткость наклонного стержня связывает перемещения по X и Y.

Исходные данные и обозначения

Узлы: (0; 0) и (3; 4) м. Ферменный элемент: A = 0,01 м², E = 200 ГПа, ρ = 7850 кг/м³. Рамный: A = 0,02 м², E = 210 ГПа, ρ = 7850 кг/м³, I = 8·10⁻⁶ м⁴. Закрепления и нагрузки не задаются: изучаются матрицы отдельного свободного элемента.

До запуска: расчёт на бумаге

1. Найдите L, c = Δx/L, s = Δy/L и осевую жёсткость стержня EA/L.
2. Вычислите три коэффициента его глобальной матрицы: K₁₁ = (EA/L)c², K₁₂ = (EA/L)cs, K₂₂ = (EA/L)s².
3. Найдите физическую массу ρAL каждой модели.

Запуск в Octave

Выполните из корня проекта. Готовая функция печатает результаты; для этой работы изменять её код не требуется.

addpath(pwd);
setup;
addpath(fullfile(pwd, "reference", "examples"));
r = example_01_element_matrices();

Работа с результатами

1. Сверьте длину и направляющие косинусы с началом вывода. В Truss local stiffness найдите осевую жёсткость, а в Truss global stiffness — три рассчитанных коэффициента. Объясните появление ненулевого K₁₂.
2. По подписям Truss DOFs и Frame DOFs подпишите строки матриц. У рамного элемента третья и шестая компоненты — повороты.
3. Рассмотрите строки 2 и 4 локальной жёсткости стержня и диагональ его матрицы масс: почему отсутствие поперечной жёсткости не означает отсутствие массы?
4. Выполните команды ниже: они находят массу при одинаковом переносе обоих узлов по X. Сверьте результат с ρAL. По формулам предскажите, как изменятся масса и осевая жёсткость при удвоении длины с сохранением направления.

tx = [1; 0; 1; 0];
fx = [1; 0; 0; 1; 0; 0];
disp(tx' * r.truss.mass * tx);
disp(fx' * r.frame.mass * fx);

Что включить в отчёт

Расчёт геометрии, таблицу трёх коэффициентов и двух масс, объяснение связи X–Y и различия между поперечной жёсткостью и массой.