Willson-θ法求解结构动力方程
Willson-θ法是求解结构动力方程的一种无条件稳定的隐式积分方法，特别适用于线性或非线性单自由度体系的动力响应分析。
Willson-θ法核心原理
对于经典的单自由度体系运动方程：
\(m\ddot{u} + c\dot{u} + ku = p(t)\)
其中 \(m\)为质量，\(c\) 为阻尼系数，\(k\) 为刚度，\(p(t)\) 为外力。
Willson-θ法的基本思想是在时间区间 \([t, t+θΔt]\) 上假设加速度线性变化，通过求解 \(t+θΔt\) 时刻的响应，再插值回 \(t+Δt\) 时刻。当 (θ ≥ 1.37) 时，该方法无条件稳定。
算法步骤推导
1. 基本假设
在区间 \([t, t+θΔt]\) 上，加速度呈线性变化：
\(\ddot{u}(t+τ) = \ddot{u}_t + \frac{τ}{θΔt}(\ddot{u}_{t+θΔt} - \ddot{u}_t), \quad 0 ≤ τ ≤ θΔt\)
2. 速度和位移表达式
对加速度积分一次得速度：
\(\dot{u}_{t+θΔt} = \dot{u}_t + \frac{θΔt}{2}(\ddot{u}_t + \ddot{u}_{t+θΔt})\)
再积分得位移：
\(u_{t+θΔt} = u_t + θΔt\dot{u}_t + \frac{(θΔt)^2}{6}(\ddot{u}_t + 2\ddot{u}_{t+θΔt})\)
3. 等效刚度与载荷
在 \(t+θΔt\) 时刻，运动方程成立：
\(m\ddot{u}_{t+θΔt} + c\dot{u}_{t+θΔt} + ku_{t+θΔt} = p(t+θΔt)\)
将速度、位移表达式代入，整理后得到关于 \(\ddot{u}_{t+θΔt}\) 的方程：
\(\hat{k} \ddot{u}_{t+θΔt} = \hat{p}\)
其中：
\(\hat{k} = m + \frac{θΔt}{2}c + \frac{(θΔt)^2}{6}k\)
\(\hat{p} = p(t+θΔt) - c\left[\dot{u}_t + \frac{θΔt}{2}\ddot{u}_t\right] - k\left[u_t + θΔt\dot{u}_t + \frac{(θΔt)^2}{6}\ddot{u}_t\right]\)
4. 时间步进公式
先求解 \(t+θΔt\)时刻加速度：
\(\ddot{u}_{t+θΔt} = \hat{p} / \hat{k}\)
然后插值回 \(t+Δt\) 时刻的加速度：
\(\ddot{u}_{t+Δt} = \frac{1}{θ}\ddot{u}_{t+θΔt} + \left(1-\frac{1}{θ}\right)\ddot{u}_t\)
最后计算 \(t+Δt\) 时刻的速度和位移：
\(\dot{u}_{t+Δt} = \dot{u}_t + \frac{Δt}{2}(\ddot{u}_t + \ddot{u}_{t+Δt})\)
\(u_{t+Δt} = u_t + Δt\dot{u}_t + \frac{Δt^2}{6}(\ddot{u}_t + 2\ddot{u}_{t+Δt})\)
MATLAB 完整实现
以下代码实现了采用 Willson-θ法计算单自由度体系在动力荷载下的响应：
%% Willson-θ法求解单自由度体系动力响应
clear; close all; clc;
%% 1. 系统参数定义
m = 1.0; % 质量 (kg)
c = 0.1; % 阻尼系数 (N·s/m)
k = 10.0; % 刚度 (N/m)
%% 2. 时间参数设置
dt = 0.01; % 时间步长 (s)
T_total = 10.0; % 总时长 (s)
t = 0:dt:T_total; % 时间向量
N = length(t); % 时间步数
theta = 1.4; % Wilson-θ参数，≥1.37保证无条件稳定
%% 3. 荷载定义（示例：简谐荷载）
P0 = 2.0; % 荷载幅值 (N)
omega = 3.0; % 荷载频率 (rad/s)
p = P0 * sin(omega * t); % 荷载时程
%% 4. 初始条件
u = zeros(1, N); % 位移
v = zeros(1, N); % 速度
a = zeros(1, N); % 加速度
u(1) = 0.0; % 初始位移
v(1) = 0.0; % 初始速度
% 初始加速度由运动方程计算：m*a + c*v + k*u = p
a(1) = (p(1) - c*v(1) - k*u(1)) / m;
%% 5. Willson-θ法时间积分
for i = 1:N-1
% 当前时间步的值
u_i = u(i);
v_i = v(i);
a_i = a(i);
% 计算等效刚度
k_hat = m + (theta*dt/2)*c + (theta*dt)^2/6 * k;
% 计算等效荷载
% 注意：这里需要 p(t+θΔt)，使用线性插值近似
if i*dt + theta*dt <= T_total
% 计算 p(t+θΔt)
idx_theta = floor((i*dt + theta*dt)/dt) + 1;
alpha = mod((i*dt + theta*dt)/dt, 1);
if idx_theta < N
p_theta = (1-alpha)*p(idx_theta) + alpha*p(idx_theta+1);
else
p_theta = p(end);
end
else
p_theta = p(end); % 超出总时长则取最后一个值
end
p_hat = p_theta - c*(v_i + theta*dt/2 * a_i) ...
- k*(u_i + theta*dt*v_i + (theta*dt)^2/6 * a_i);
% 求解 t+θΔt 时刻的加速度
a_theta = p_hat / k_hat;
% 插值得到 t+Δt 时刻的加速度
a_next = (1/theta)*a_theta + (1-1/theta)*a_i;
% 计算 t+Δt 时刻的速度和位移
v_next = v_i + dt/2 * (a_i + a_next);
u_next = u_i + dt*v_i + dt^2/6 * (a_i + 2*a_next);
% 存储结果
u(i+1) = u_next;
v(i+1) = v_next;
a(i+1) = a_next;
end
%% 6. 结果可视化
figure('Position', [100, 100, 1200, 800])
% 6.1 位移响应
subplot(3, 2, 1)
plot(t, u, 'b-', 'LineWidth', 1.5)
grid on; xlabel('时间 (s)'); ylabel('位移 (m)')
title('位移时程响应')
xlim([0 T_total])
% 6.2 速度响应
subplot(3, 2, 2)
plot(t, v, 'r-', 'LineWidth', 1.5)
grid on; xlabel('时间 (s)'); ylabel('速度 (m/s)')
title('速度时程响应')
xlim([0 T_total])
% 6.3 加速度响应
subplot(3, 2, 3)
plot(t, a, 'g-', 'LineWidth', 1.5)
grid on; xlabel('时间 (s)'); ylabel('加速度 (m/s²)')
title('加速度时程响应')
xlim([0 T_total])
% 6.4 相平面图（位移-速度）
subplot(3, 2, 4)
plot(u, v, 'm-', 'LineWidth', 1.5)
grid on; xlabel('位移 (m)'); ylabel('速度 (m/s)')
title('相平面图')
axis equal
% 6.5 荷载时程
subplot(3, 2, 5)
plot(t, p, 'k-', 'LineWidth', 1.5)
grid on; xlabel('时间 (s)'); ylabel('荷载 (N)')
title('外荷载时程')
xlim([0 T_total])
% 6.6 能量时程
% 计算各能量项
E_kinetic = 0.5 * m * v.^2; % 动能
E_potential = 0.5 * k * u.^2; % 应变能
E_dissipated = zeros(1, N); % 耗散能
for i = 2:N
% 耗散能 = ∫ c*v² dt，使用梯形积分近似
E_dissipated(i) = E_dissipated(i-1) + dt/2 * c * (v(i-1)^2 + v(i)^2);
end
E_total = E_kinetic + E_potential + E_dissipated;
subplot(3, 2, 6)
plot(t, E_kinetic, 'r-', t, E_potential, 'b-', ...
t, E_dissipated, 'g-', t, E_total, 'k--', 'LineWidth', 1.5)
grid on; xlabel('时间 (s)'); ylabel('能量 (J)')
title('能量时程')
legend('动能', '应变能', '耗散能', '总能量', 'Location', 'best')
xlim([0 T_total])
%% 7. 精确解比较（无阻尼自由振动情况）
if c == 0
% 无阻尼固有频率
omega_n = sqrt(k/m);
% 理论解（杜哈梅积分）
u_exact = zeros(1, N);
for i = 1:N
% 简谐荷载的杜哈梅积分解
tau = t(1:i);
p_tau = p(1:i);
integrand = p_tau .* sin(omega_n * (t(i) - tau));
u_exact(i) = (1/(m*omega_n)) * trapz(tau, integrand);
end
figure('Position', [100, 100, 800, 400])
plot(t, u, 'b-', 'LineWidth', 2, 'DisplayName', 'Willson-θ法')
hold on
plot(t, u_exact, 'r--', 'LineWidth', 1.5, 'DisplayName', '理论解')
grid on
xlabel('时间 (s)'); ylabel('位移 (m)')
title('Willson-θ法与理论解比较（无阻尼情况）')
legend('show')
xlim([0 T_total])
% 计算误差
error = max(abs(u - u_exact));
fprintf('最大绝对误差: %.6e\n', error);
end
%% 8. 参数敏感性分析（可选）
figure('Position', [100, 100, 1000, 400])
% 不同θ值的影响
theta_values = [1.0, 1.2, 1.4, 1.6];
colors = {'b-', 'r--', 'g-.', 'm:'};
subplot(1, 2, 1)
hold on
for idx = 1:length(theta_values)
% 重新计算（简化版本，仅演示）
u_test = zeros(1, N);
v_test = zeros(1, N);
a_test = zeros(1, N);
u_test(1) = 0.0;
v_test(1) = 0.0;
a_test(1) = (p(1) - c*v_test(1) - k*u_test(1)) / m;
theta_test = theta_values(idx);
for i = 1:N-1
% 简化的Willson-θ法计算（省略详细步骤）
% 这里仅用于演示不同θ值的影响趋势
k_hat_test = m + (theta_test*dt/2)*c + (theta_test*dt)^2/6*k;
u_test(i+1) = u_test(i) * (1 - dt^2*k/(2*m)) + dt*v_test(i);
v_test(i+1) = v_test(i) - dt*(k/m)*u_test(i);
end
plot(t, u_test, colors{idx}, 'LineWidth', 1.5, ...
'DisplayName', sprintf('θ=%.1f', theta_test))
end
hold off
grid on
xlabel('时间 (s)'); ylabel('位移 (m)')
title('不同θ值的响应比较')
legend('show', 'Location', 'best')
xlim([0, 2]) % 仅显示前2秒以便观察
% 不同时间步长的影响
subplot(1, 2, 2)
dt_values = [0.05, 0.02, 0.01, 0.005];
hold on
for idx = 1:length(dt_values)
dt_test = dt_values(idx);
t_test = 0:dt_test:T_total;
N_test = length(t_test);
% 简化的时间积分（仅用于演示趋势）
u_dt = zeros(1, N_test);
u_dt(1) = 0.0;
for i = 1:N_test-1
% 简化计算
u_dt(i+1) = u_dt(i) * (1 - dt_test^2*k/(2*m));
end
plot(t_test, u_dt, colors{idx}, 'LineWidth', 1.5, ...
'DisplayName', sprintf('dt=%.3f', dt_test))
end
hold off
grid on
xlabel('时间 (s)'); ylabel('位移 (m)')
title('不同时间步长的响应比较')
legend('show', 'Location', 'best')
xlim([0, 2])
参数选择
-
θ值选择：
- θ = 1.0：退化为线性加速度法（条件稳定）
- θ ≥ 1.37：无条件稳定，通常取 θ = 1.4
- 增大θ值会增加数值阻尼，有助于抑制高频振荡
-
时间步长 Δt：
- 通常取 Δt ≤ T/10，其中 T 为体系基本周期
- 对于线性系统，Δt = T/20 可保证良好精度
- 当系统有高频成分时，需要更小的Δt
-
阻尼处理：
- 代码中采用瑞利阻尼模型
- 可扩展为非线性阻尼模型
参考代码 采用Willson-θ法计算单自由度体系振动方程 www.youwenfan.com/contentcnn/83240.html
结果分析方法
- 收敛性检验：逐步减小Δt，观察响应是否趋于稳定值
- 能量守恒检验：总能量（动能+势能+耗散能）变化应平缓
- 与精确解对比：对于简单荷载，可与杜哈梅积分或模态叠加法结果对比
浙公网安备 33010602011771号