MATLAB实现Bouc-Wen模型动力响应计算
Bouc-Wen模型是一种广泛使用的非线性滞回模型，用于描述材料在循环加载下的非线性滞回行为。
1. Bouc-Wen模型基本原理
Bouc-Wen模型的基本方程如下：
m*x'' + c*x' + F(x, z) = f(t)
z' = A*x' - β*|x'|*|z|^(n-1)*z - γ*x'*|z|^n
F(x, z) = α*k*x + (1-α)*k*z
其中：
- m: 质量
- c: 阻尼系数
- k: 刚度系数
- α: 屈服后刚度与初始刚度之比
- A, β, γ, n: Bouc-Wen模型参数
- x: 位移
- z: 滞回变量
- f(t): 外部激励
2. MATLAB实现代码
function bouc_wen_simulation()
% Bouc-Wen模型参数
params.m = 1; % 质量 (kg)
params.c = 0.2; % 阻尼系数 (N·s/m)
params.k = 100; % 初始刚度 (N/m)
params.alpha = 0.1; % 屈服后刚度与初始刚度之比
params.A = 1; % Bouc-Wen参数
params.beta = 0.5; % Bouc-Wen参数
params.gamma = 0.5; % Bouc-Wen参数
params.n = 1; % Bouc-Wen参数
% 仿真参数
tspan = [0, 10]; % 时间范围 (s)
dt = 0.01; % 时间步长 (s)
t = tspan(1):dt:tspan(2);
% 外部激励 (正弦波)
f_amp = 10; % 激励幅值 (N)
f_freq = 1; % 激励频率 (Hz)
f_t = f_amp * sin(2*pi*f_freq*t);
% 初始条件
x0 = 0; % 初始位移 (m)
v0 = 0; % 初始速度 (m/s)
z0 = 0; % 初始滞回变量
% 使用ode45求解微分方程
options = odeset('RelTol', 1e-6, 'AbsTol', 1e-9);
[T, Y] = ode45(@(t, y) bouc_wen_ode(t, y, f_t, tspan, params), tspan, [x0; v0; z0], options);
% 插值以获得均匀时间步长的结果
t_uniform = linspace(tspan(1), tspan(2), length(t));
x = interp1(T, Y(:,1), t_uniform);
v = interp1(T, Y(:,2), t_uniform);
z = interp1(T, Y(:,3), t_uniform);
% 计算恢复力
F = zeros(size(x));
for i = 1:length(x)
F(i) = params.alpha*params.k*x(i) + (1-params.alpha)*params.k*z(i);
end
% 绘制结果
plot_results(t_uniform, x, v, z, F, f_t, params);
end
function dydt = bouc_wen_ode(t, y, f_t, tspan, params)
% 从y中提取状态变量
x = y(1); % 位移
v = y(2); % 速度
z = y(3); % 滞回变量
% 插值获取当前时间的激励
t_vec = linspace(tspan(1), tspan(2), length(f_t));
f = interp1(t_vec, f_t, t);
% Bouc-Wen模型方程
dzdt = params.A*v - params.beta*abs(v)*abs(z)^(params.n-1)*z - params.gamma*v*abs(z)^params.n;
% 恢复力
F = params.alpha*params.k*x + (1-params.alpha)*params.k*z;
% 运动方程
dxdt = v;
dvdt = (f - params.c*v - F) / params.m;
% 返回导数
dydt = [dxdt; dvdt; dzdt];
end
function plot_results(t, x, v, z, F, f_t, params)
% 创建图形
figure('Position', [100, 100, 1200, 800]);
% 位移时程
subplot(3, 2, 1);
plot(t, x, 'b', 'LineWidth', 1.5);
xlabel('时间 (s)');
ylabel('位移 (m)');
title('位移时程');
grid on;
% 速度时程
subplot(3, 2, 2);
plot(t, v, 'r', 'LineWidth', 1.5);
xlabel('时间 (s)');
ylabel('速度 (m/s)');
title('速度时程');
grid on;
% 滞回变量时程
subplot(3, 2, 3);
plot(t, z, 'g', 'LineWidth', 1.5);
xlabel('时间 (s)');
ylabel('滞回变量 z');
title('滞回变量时程');
grid on;
% 恢复力时程
subplot(3, 2, 4);
plot(t, F, 'm', 'LineWidth', 1.5);
xlabel('时间 (s)');
ylabel('恢复力 (N)');
title('恢复力时程');
grid on;
% 滞回曲线 (力-位移关系)
subplot(3, 2, 5);
plot(x, F, 'k', 'LineWidth', 1.5);
xlabel('位移 (m)');
ylabel('恢复力 (N)');
title('滞回曲线');
grid on;
% 激励时程
subplot(3, 2, 6);
plot(t, f_t, 'c', 'LineWidth', 1.5);
xlabel('时间 (s)');
ylabel('激励力 (N)');
title('外部激励');
grid on;
% 添加参数信息
annotation('textbox', [0.15, 0.01, 0.7, 0.05], 'String', ...
sprintf('参数: m=%.1f, c=%.1f, k=%.1f, α=%.1f, A=%.1f, β=%.1f, γ=%.1f, n=%.1f', ...
params.m, params.c, params.k, params.alpha, params.A, params.beta, params.gamma, params.n), ...
'FitBoxToText', 'on', 'EdgeColor', 'none', 'FontSize', 10);
end
3. 参数敏感性分析
为了研究不同参数对Bouc-Wen模型响应的影响，我们可以编写一个参数敏感性分析函数：
function parameter_sensitivity_analysis()
% 基准参数
base_params.m = 1;
base_params.c = 0.2;
base_params.k = 100;
base_params.alpha = 0.1;
base_params.A = 1;
base_params.beta = 0.5;
base_params.gamma = 0.5;
base_params.n = 1;
% 仿真参数
tspan = [0, 10];
dt = 0.01;
t = tspan(1):dt:tspan(2);
f_amp = 10;
f_freq = 1;
f_t = f_amp * sin(2*pi*f_freq*t);
% 初始条件
x0 = 0;
v0 = 0;
z0 = 0;
% 测试不同参数值
alpha_values = [0.05, 0.1, 0.2];
beta_values = [0.3, 0.5, 0.7];
gamma_values = [0.3, 0.5, 0.7];
n_values = [1, 2, 3];
% 创建图形
figure('Position', [100, 100, 1200, 1000]);
% α参数敏感性分析
subplot(2, 2, 1);
hold on;
for i = 1:length(alpha_values)
params = base_params;
params.alpha = alpha_values(i);
[T, Y] = ode45(@(t, y) bouc_wen_ode(t, y, f_t, tspan, params), tspan, [x0; v0; z0]);
x = interp1(T, Y(:,1), t);
z = interp1(T, Y(:,3), t);
F = params.alpha*params.k*x + (1-params.alpha)*params.k*z;
plot(x, F, 'LineWidth', 1.5, 'DisplayName', sprintf('α=%.2f', params.alpha));
end
xlabel('位移 (m)');
ylabel('恢复力 (N)');
title('α参数敏感性分析');
legend('show');
grid on;
hold off;
% β参数敏感性分析
subplot(2, 2, 2);
hold on;
for i = 1:length(beta_values)
params = base_params;
params.beta = beta_values(i);
[T, Y] = ode45(@(t, y) bouc_wen_ode(t, y, f_t, tspan, params), tspan, [x0; v0; z0]);
x = interp1(T, Y(:,1), t);
z = interp1(T, Y(:,3), t);
F = params.alpha*params.k*x + (1-params.alpha)*params.k*z;
plot(x, F, 'LineWidth', 1.5, 'DisplayName', sprintf('β=%.1f', params.beta));
end
xlabel('位移 (m)');
ylabel('恢复力 (N)');
title('β参数敏感性分析');
legend('show');
grid on;
hold off;
% γ参数敏感性分析
subplot(2, 2, 3);
hold on;
for i = 1:length(gamma_values)
params = base_params;
params.gamma = gamma_values(i);
[T, Y] = ode45(@(t, y) bouc_wen_ode(t, y, f_t, tspan, params), tspan, [x0; v0; z0]);
x = interp1(T, Y(:,1), t);
z = interp1(T, Y(:,3), t);
F = params.alpha*params.k*x + (1-params.alpha)*params.k*z;
plot(x, F, 'LineWidth', 1.5, 'DisplayName', sprintf('γ=%.1f', params.gamma));
end
xlabel('位移 (m)');
ylabel('恢复力 (N)');
title('γ参数敏感性分析');
legend('show');
grid on;
hold off;
% n参数敏感性分析
subplot(2, 2, 4);
hold on;
for i = 1:length(n_values)
params = base_params;
params.n = n_values(i);
[T, Y] = ode45(@(t, y) bouc_wen_ode(t, y, f_t, tspan, params), tspan, [x0; v0; z0]);
x = interp1(T, Y(:,1), t);
z = interp1(T, Y(:,3), t);
F = params.alpha*params.k*x + (1-params.alpha)*params.k*z;
plot(x, F, 'LineWidth', 1.5, 'DisplayName', sprintf('n=%d', params.n));
end
xlabel('位移 (m)');
ylabel('恢复力 (N)');
title('n参数敏感性分析');
legend('show');
grid on;
hold off;
end
4. 使用不同激励函数
Bouc-Wen模型可以响应不同类型的激励。以下是如何使用不同激励函数的示例：
function different_excitations()
% Bouc-Wen模型参数
params.m = 1;
params.c = 0.2;
params.k = 100;
params.alpha = 0.1;
params.A = 1;
params.beta = 0.5;
params.gamma = 0.5;
params.n = 1;
% 仿真参数
tspan = [0, 10];
dt = 0.01;
t = tspan(1):dt:tspan(2);
% 初始条件
x0 = 0;
v0 = 0;
z0 = 0;
% 创建不同激励
% 1. 正弦波
f_sine = 10 * sin(2*pi*1*t);
% 2. 白噪声
rng(0); % 设置随机种子以确保可重复性
f_noise = 5 * randn(size(t));
% 3. 地震波 (简化模拟)
f_earthquake = zeros(size(t));
earthquake_start = find(t >= 2, 1);
earthquake_end = find(t >= 4, 1);
f_earthquake(earthquake_start:earthquake_end) = ...
15 * sin(2*pi*2*(t(earthquake_start:earthquake_end)-2)) .* ...
exp(-0.5*(t(earthquake_start:earthquake_end)-2));
% 计算响应
excitations = {f_sine, f_noise, f_earthquake};
names = {'正弦激励', '白噪声激励', '地震激励'};
% 绘制结果
figure('Position', [100, 100, 1200, 800]);
for i = 1:3
f_t = excitations{i};
% 求解微分方程
[T, Y] = ode45(@(t, y) bouc_wen_ode(t, y, f_t, tspan, params), tspan, [x0; v0; z0]);
% 插值以获得均匀时间步长的结果
t_uniform = linspace(tspan(1), tspan(2), length(t));
x = interp1(T, Y(:,1), t_uniform);
z = interp1(T, Y(:,3), t_uniform);
F = params.alpha*params.k*x + (1-params.alpha)*params.k*z;
% 绘制滞回曲线
subplot(2, 3, i);
plot(x, F, 'b', 'LineWidth', 1.5);
xlabel('位移 (m)');
ylabel('恢复力 (N)');
title(sprintf('%s - 滞回曲线', names{i}));
grid on;
% 绘制激励时程
subplot(2, 3, i+3);
plot(t, f_t, 'r', 'LineWidth', 1.5);
xlabel('时间 (s)');
ylabel('激励力 (N)');
title(sprintf('%s - 激励时程', names{i}));
grid on;
end
end
5. 如何使用这些函数
要运行上述代码，只需在MATLAB命令窗口中调用相应的函数：
% 运行基本Bouc-Wen模型仿真
bouc_wen_simulation();
% 运行参数敏感性分析
parameter_sensitivity_analysis();
% 运行不同激励下的响应分析
different_excitations();
参考代码 用于计算bouc-wen模型的动力响应 www.3dddown.com/cna/50851.html
6. 事项
-
参数选择：Bouc-Wen模型的参数需要根据具体材料或结构的实验数据来确定。不同参数组合会产生不同的滞回行为。
-
数值稳定性：对于某些参数组合，特别是当n值较大时，模型可能会变得数值不稳定。可能需要减小时间步长或使用更高级的数值积分方法。
-
初始条件：对于强非线性系统，响应可能对初始条件敏感。在实际应用中，可能需要考虑不同的初始条件。
-
模型扩展：标准Bouc-Wen模型可以扩展为更复杂的版本，如考虑刚度退化、强度退化等效应的扩展Bouc-Wen模型。
浙公网安备 33010602011771号