Два элемента, обе матрицы, три формы

1. Расчётная схема

Два призматических рамных элемента соединяются в узле 2. Узлы 1 и 3 полностью закреплены, поэтому из девяти глобальных степеней свободы остаются только три компоненты среднего узла: u_x2, u_y2 и θ₂.

Точные координаты наклонного элемента равны (100/√2, 100/√2). Поэтому входной файл использует 70.7106781186548, а не сокращённое 70.7: последнее изменило бы длину первого элемента на 0.0151%.

2. Исполняемая модель

Вся схема помещается в короткий входной файл. Нагрузка оставлена как часть исходного примера, но метод RunModal её не использует.

nodes
3
0,0,0
70.7106781186548,70.7106781186548,0
170.710678118655,70.7106781186548,0
elems_113
2
113,1,2,6,10000000,0.7,100
113,2,3,6,10000000,0.7,100
bcfix
2
1,1,0,0,0
1,3,0,0,0
bcforce_stat
1
10,2,0,-10,0

Для рамного узла карта СС имеет вид [ux, uy, thetaZ]. Поэтому свободные СС узла 2 получают глобальные номера [4 5 6], а закреплённые СС — [1 2 3 7 8 9].

3. Локальные матрицы элемента

Оба элемента имеют одинаковые свойства и длину, поэтому их локальные матрицы совпадают. Численные коэффициенты жёсткости:

Локальная матрица жёсткости, записанная в порядке [u1 v1 θ1 u2 v2 θ2]:

k(local) = 10^6 ×
[ 0.600   0       0      -0.600   0       0
  0       0.012   0.600   0      -0.012   0.600
  0       0.600  40.000   0      -0.600  20.000
 -0.600   0       0       0.600   0       0
  0      -0.012  -0.600   0       0.012  -0.600
  0       0.600  20.000   0      -0.600  40.000 ]

Для согласованной массы множитель перед безразмерным шаблоном особенно удобен:

ρAL/420 = 0.7·6·100/420 = 1

m(local) =
[ 140      0       0      70      0       0
    0    156    2200       0     54   -1300
    0   2200   40000       0   1300  -30000
   70      0       0     140      0       0
    0     54    1300       0    156   -2200
    0  -1300  -30000       0  -2200   40000 ]

4. Поворот в глобальные оси

Для первого элемента c=s=1/√2. Книга выносит общий множитель 1/√2; реализация записывает те же направляющие косинусы непосредственно:

T1 = 1/sqrt(2) ×
[ 1   1   0   0   0   0
 -1   1   0   0   0   0
  0   0  sqrt(2)  0   0   0
  0   0   0   1   1   0
  0   0   0  -1   1   0
  0   0   0   0   0  sqrt(2) ]

T2 = I

Преобразование применяется одинаково к жёсткости и массе:

k(global) = T^T k(local) T
m(global) = T^T m(local) T

Именно здесь наклонный элемент создаёт связи между глобальными компонентами x, y и поворотом. Горизонтальный элемент преобразования не требует.

5. Сборка в трёх свободных координатах

После закрепления наружных узлов достаточно сложить свободно-свободные блоки двух элементов. Это тот же прямой метод, который для большой модели складывает вклады по глобальной карте СС.

Вклады элемента 1, наклон 45°

K1(ff) =
[ 306000   294000    424264.0687
  294000   306000   -424264.0687
  424264.0687 -424264.0687 40000000 ]

M1(ff) =
[ 148    -8     1555.6349
   -8   148    -1555.6349
 1555.6349 -1555.6349 40000 ]

Вклады элемента 2, горизонталь

K2(ff) =
[ 600000       0         0
       0   12000    600000
       0  600000  40000000 ]

M2(ff) =
[ 140    0      0
    0  156   2200
    0 2200  40000 ]

Редуцированные матрицы системы

Kff = K1(ff) + K2(ff)
    = 10^6 ×
[ 0.906000  0.294000  0.424264
  0.294000  0.318000  0.175736
  0.424264  0.175736 80.000000 ]

Mff = M1(ff) + M2(ff)
    =
[ 288        -8      1555.635
   -8       304       644.365
 1555.635   644.365 80000 ]

После округления 0.424264 → 0.424, 0.175736 → 0.176, 1555.635 → 1556 и 644.365 → 644 получаются ровно матрицы, напечатанные в книге.

6. Собственные значения и частоты

Ненулевое решение существует, когда:

det(K_ff − ω²M_ff) = 0

В численном виде решается обобщённая задача eig(Kff, Mff). Её собственные значения равны λ=ω², после чего вычисляются ω=√λ и f=ω/(2π).

7. Массово-нормированные формы

Для каждой формы выбирается масштаб, при котором φ_i^TMφ_i=1. Столбцы модальной матрицы:

Phi =
[  0.02182975   0.00497889   0.05830737
  -0.05270167   0.00206232   0.02415170
   0            0.00340929  -0.00162915 ]

Напечатанные в книге значения [0.0218, -0.0527, 0], [0.00498, 0.00206, 0.00341] и [0.0583, 0.0241, -0.0016] являются округлением этих столбцов.

Две обязательные численные проверки:

Φ^TM_ffΦ = I

r_i = K_ffφ_i − λ_iM_ffφ_i ≈ 0

Знак собственного вектора произволен: φ и −φ описывают одну форму. Поэтому тест сначала нормирует каждый столбец по массе, затем согласует его знак с эталоном и лишь после этого сравнивает компоненты.

8. Воспроизведение в GNU Octave

addpath(pwd);
setup;

problem = StructFEProblem( ...
    fullfile('examples', 'cases', ...
        'CaseFig11.7p363 MarioPaz.txt'), ...
    struct('verbose', false, 'plotting', false));

model = problem.GetAnalysisModel();
result = problem.RunModal();

Kff = full(model.stiffness(result.freeDOFs, result.freeDOFs))
Mff = full(model.mass(result.freeDOFs, result.freeDOFs))
lambda = result.eigenvalues
omega = result.angularFrequenciesRadPerSec
frequencyHz = result.frequenciesHz
Phi = result.modeShapes(result.freeDOFs, :)

Автоматическая проверка находится в tests/test_mario_paz_modal.m и запускается вместе со всем набором:

addpath(pwd); setup; run_octave_tests;

9. Практические эксперименты