Метод средних ускорений

Уравнение движения и редукция

На свободных СС решается система:

M_ff a(t) + K_ff u(t) = F_f(t)

Это линейный недемпфированный частный случай более общего уравнения с матрицей демпфирования C и нелинейной восстанавливающей силой. В МКЭ‑Ф C = 0, нелинейного слагаемого нет, а закреплённые компоненты полных историй остаются точно нулевыми.

Единая временная сетка

При шаге Δt и запрошенной длительности T число шагов равно fix(T/dt). Все возвращаемые истории имеют N + 1 столбец:

t₀ = 0,   t₁ = Δt,   …,   t_N = NΔt

Если T/Δt нецелое, фактическая длительность result.duration меньше запрошенной result.requestedDuration. Номер столбца в Octave на единицу больше номера шага: столбец 1 соответствует t = 0.

Начальное состояние

Пользователь может задать полные векторы u₀ и v₀. Их закреплённые компоненты обязаны быть нулевыми. Начальное ускорение не задаётся независимо, а вычисляется из равновесия:

a₀ = M_ff⁻¹[F_f(0) − K_ffu₀]

Так три начальные истории начинают расчёт из одного физически согласованного состояния.

Аппроксимации Ньюмарка

Перемещение и скорость на следующем шаге связываются с ускорениями:

u_n+1 = u_n + Δt v_n + Δt²[(1/2 − β)a_n + βa_n+1]

v_n+1 = v_n + Δt[(1 − γ)a_n + γa_n+1]

Выбор β = 1/4, γ = 1/2 задаёт метод средних ускорений. Для линейной системы он имеет второй порядок точности и безусловно устойчив, но устойчивость не означает достаточную точность при слишком крупном Δt.

Эффективная жёсткость

В коде используются коэффициенты:

a0 = 1/(beta*dt^2);
a2 = 1/(beta*dt);
a3 = 1/(2*beta) - 1;

На каждом шаге решается одна система:

K_eff u_n+1 = F_n+1 + M(a₀u_n + a₂v_n + a₃a_n)

K_eff = K_ff + a₀M_ff

Матрица K_eff постоянна, поэтому solveTransient один раз выполняет разложение Холецкого и повторно использует его на всех шагах. После решения обновляются ускорение и скорость.

Импульс на дискретной сетке

Нагрузка типа 12 ненулевая только в столбце t = Δt. Для прямоугольного одношагового импульса:

I_discrete = F₀Δt

Поэтому при сравнении сеток с разным Δt следует сохранять не амплитуду, а импульс: F₀ = I/Δt. Иначе уточнение шага одновременно меняет физическое возбуждение и численную погрешность.

Сходимость и энергия

Для свободного недемпфированного осциллятора с ω = √(k/m):

u(t) = u₀cos(ωt) + [v₀/ω]sin(ωt)

Уменьшение Δt вдвое должно примерно в четыре раза уменьшать ошибку метода второго порядка. Это надёжнее одной проверки при одном шаге.

Полная механическая энергия свободной системы:

E(t) = 1/2 v^TMv + 1/2 u^TKu

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

Односторонний спектр

computeSingleSidedSpectrum применяет БПФ к каждой строке перемещений и возвращает неотрицательную половину частот. При N_s отсчётах:

Δf = 1/(N_sΔt)

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

Как запустить расчёт и посмотреть результат

В Octave из корневого каталога проекта выполните:

addpath(pwd);
setup;

problem = StructFEProblem(fullfile( ...
    'examples', 'cases', 'CaseBeam.txt'));
model = problem.GetAnalysisModel();
result = problem.RunTransient(1e-4, 1e-2, 3, 2);

options = struct(...
    'title', 'Переходный отклик', ...
    'lengthUnit', 'm', 'forceUnit', 'N', ...
    'momentUnit', 'N*m', 'timeUnit', 's', ...
    'transientFields', {{'displacements', 'velocities', ...
                         'accelerations', 'loadHistory', ...
                         'reactions', 'spectrum'}}, ...
    'timeStride', 1, 'prettyPrint', true);
exportPostprocessorData(model, result, ...
    fullfile('output', 'transient-history.json'), options);

Затем откройте postprocessor/index.html двойным щелчком в Chrome и перетащите на страницу output/transient-history.json или выберите его через Open JSON. Готовый пример находится в examples/results/transient-history.json. Локальный сервер не нужен.

Лабораторная 06. Свободные колебания: точность метода Ньюмарка

Цель — сравнить численное движение с аналитическим и определить по трём сеткам порядок сходимости.

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

m = 2 кг, k = 18 Н/м, u₀ = 0,02 м, v₀ = −0,01 м/с. Сила и демпфирование равны нулю. Готовый пример рассчитывает три периода с 50, 100 и 200 шагами на период; β = 1/4, γ = 1/2.

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

1. Найдите ω = √(k/m), f = ω/(2π), T = 1/f и начальное ускорение a₀ = −ku₀/m.
2. Запишите u(t) = u₀cos(ωt) + (v₀/ω)sin(ωt). Найдите начальную энергию E₀ = mv₀²/2 + ku₀²/2.

Запуск в Octave

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

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

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

1. Сверьте Natural frequency, period и Initial acceleration с расчётом.
2. Из Steps per period and maximum displacement error выпишите три ошибки. Рассчитайте отношения ошибки на 50 шагах к ошибке на 100 и на 100 к 200. Почему отношение около четырёх указывает на второй порядок?
3. Готовыми командами ниже постройте точное и численное движение на самой мелкой сетке, а затем их разность. Отдельный график разности позволяет увидеть погрешность при визуально совпадающих кривых.
4. Сопоставьте Relative energy drift с найденной погрешностью движения. Объясните, почему сохранение энергии ещё не означает точного совпадения фазы.
5. Сверьте Spectrum peak с собственной частотой: расхождение должно укладываться в напечатанную bin width. По формуле Δf = 1/(NΔt) предскажите влияние сокращения длительности записи при прежнем шаге.

После запуска выполните в Octave с доступным графическим окном:

figure;
plot(r.time, r.displacements(1,:), "-", ...
     r.time, r.exactDisplacement, "--");
xlabel("t [s]"); ylabel("u [m]");
legend("Newmark", "Exact"); grid on;
figure;
plot(r.time, r.displacements(1,:) - r.exactDisplacement);
xlabel("t [s]"); ylabel("u - u_exact [m]"); grid on;

При безоконном запуске вместо графиков выпишите несколько пар значений и их разность:

idx = round(linspace(1, numel(r.time), 13));
disp([r.time(idx).', r.displacements(1,idx).', ...
      r.exactDisplacement(idx).', ...
      (r.displacements(1,idx) - r.exactDisplacement(idx)).']);

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

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

MPC

Newmark интегрирует независимые координаты с Kr, Mr и TᵀF(t). Перемещения, скорости и ускорения восстанавливаются через T. Начальные перемещения и скорости должны удовлетворять SPC и MPC; несовместные значения отклоняются. Динамические силы связей восстанавливаются из r = Ma + Ku − F.

Подробное описание MPC · JSON v3 · CaseMPCFrame.txt · CaseMPCTruss.txt