一、背景
研究振动体系对于周期荷载的结构响应状态,可借助傅里叶变换,将周期荷载
转换为三角波的线性叠加后,运用基本的结构动力学知识即可解出体系位移情况。
本文重点内容:
- 分段函数的傅里叶展开
- 计算阻尼对振动的影响
二、理论推导
一般谐振荷载下的振动方程:
对周期荷载进行傅里叶展开:
其中的傅里叶系数:
1.无阻尼情况
解出的通解:
其中
2.考虑阻尼情况
解出的通解:
其中
三、代码实现
3.1物理模型
3.2实现分段函数的积分
此处默认外力的分段函数为两段,integrals文件的完整代码
function [ output ] = integrals( ... % 此处默认外力的分段函数为两段
t, ... % 积分变量
fx_1, ... % 第一段函数表达式
leftInterval_1, ... % 左区间
rightInterval_1, ... % 右区间
fx_2, ... % 第二段函数表达式
leftInterval_2, ... % 左区间
rightInterval_2 ... % 右区间
)
output=int(fx_1,t,leftInterval_1,rightInterval_1) + ... % 计算第一段函数的积分
int(fx_2,t,leftInterval_2, rightInterval_2); % 计算第二段函数的积分
end
3.3傅里叶展开外力函数、无阻尼时的位移函数、有阻尼时的位移函数
method_fourier_unfolds文件的完整代码
function [ ...
external_force, ... % 经傅里叶展开后的外力函数
displacement_nodump, ... % 无阻尼时的位移函数
displacement_dump ... % 存在阻尼时的位移函数
] = method_fourier_unfolds( ...
t1, ... % 前半周期时长
t2, ... % 后半周期时长
exact_external_force_1, ...
exact_external_force_2, ...
expand_num, ... % 傅里叶展开项数
omega, ... % 体系固有频率
k, ... % 结构刚度系数
kesai ... % 阻尼比
)
T = t1 + t2; % 荷载总周期
theta=2*pi/T; % 外荷载频率
% 分段函数进行傅里叶展开的核心代码
syms t n;
a0=1/T * integrals('t',exact_external_force_1,0,t1,exact_external_force_2,t1,T);
an=2/T * integrals('t',exact_external_force_1*cos(n*theta*t),0,t1,exact_external_force_2*cos(n*theta*t),t1,T);
bn=2/T * integrals('t',exact_external_force_1*sin(n*theta*t),0,t1,exact_external_force_2*sin(n*theta*t),t1,T);
% 初始化必要的值
displacement_nodump = 0;
displacement_dump = 0;
beta = 1:expand_num;
miu = 1:expand_num;
external_force = 0;
% 以传入的傅里叶展开项数进行计算,最后将结果叠加即是傅里叶展开的近似
for n=1:expand_num
% 以下运算没有复杂的逻辑,也未涉及复杂的循环迭代,只需要按数学表达式写出即可
beta(n) = n*theta/omega;
miu(n) = abs(1./(1-beta(n).^2));
epsilon = atan(2*kesai*beta(n)/(1-beta(n).^2));
miu_dump = 1./((1-beta(n).^2).^2+(2*kesai*beta(n)).^2).^0.5;
displacement_nodump=displacement_nodump + miu(n) .* vpa(eval(an*cos(n*theta.*t)+bn*sin(n*theta.*t)));
if (nargin > 7 ) % 若输入值大于7,说明需要同时计算有阻尼和无阻尼情况,否则只计算无阻尼情况
displacement_dump = displacement_dump + miu_dump .* vpa(eval(an*cos(n*theta.*t-epsilon)+bn*sin(n*theta.*t-epsilon)));
end
external_force = external_force+vpa(eval(an)*cos(n*theta.*t)+eval(bn)*sin(n*theta.*t));
end
external_force = external_force+a0; % 返回作用力P的表达式
displacement_nodump = (displacement_nodump+double(a0))/k; % 返回无阻尼时的位移表达式
displacement_dump = (displacement_dump+double(a0))/k; % 返回有阻尼时的位移表达式
end
3.4main文件中最终实现(完整代码)
clear all
%% 设置参数
m = 2e4; % 小球质量
EI = 1e6; % 刚度
L = 5; % 杆长
k = 3*EI/(L^3); % 计算体系刚度系数
w = (k/m)^0.5; % 计算体系固有频率
p_max = 1000; % 最大外力值
t1 = 1; % 前半周期
t2 = 1; % 后半周期
% 定义符合变量t
syms t
% 方波荷载
force_1 = 0*t+p_max;
force_2 = 0*t;
% 进行傅里叶展开
[external_force, displacement_nodump, displacement_dump] = ...
method_fourier_unfolds(t1, t2, force_1, force_2, 50, w, k, 0.6);
len = 200;
y1 = 1:len;
y2 = 1:len;
p = 1:len;
time = 1:len;
t = 0;
for i=1:len
t = t + 0.05;
time(i) = t;
y1(i) = eval(displacement_nodump);
y2(i) = eval(displacement_dump);
p(i) = eval(external_force);
end
subplot(2, 2, 1);
plot(time, y1);
title('无阻尼时位移与时间曲线');
subplot(2, 2, 2);
plot(time, y2);
title('有阻尼时位移与时间曲线');
subplot(2, 2, [3,4]);
plot(time, p);
title('傅里叶展开的外力函数图像');
结果
四、补充说明
因为振动方程中考虑了结构的众多物理力学性质,如质量,弹性模量,刚度、固有频率等,而本文借助matlab实现时,未对公式做简化,均采用变量形式进行运算,所以可以在此源码基础上,进一步探讨不同物理量对结构振动的影响。(以下为阻尼比对结构振动的影响)
如有任何疑问,请在评论区留言,
源文件github地址:https://github.com/darlingxyz/method_fourier_expansion_in_structural_mechanics