【Matlab】结构在傅里叶展开下的周期荷载响应——文末附源码

最新推荐文章于 2023-04-27 19:06:52 发布

4399资深专家

最新推荐文章于 2023-04-27 19:06:52 发布

阅读量688

点赞数 1

分类专栏：数学物理编程文章标签： matlab 算法

本文链接：https://blog.csdn.net/qq_50920297/article/details/124236017

版权

数学物理编程专栏收录该内容

11 篇文章 15 订阅

订阅专栏

一、背景

研究振动体系对于周期荷载的结构响应状态，可借助傅里叶变换，将周期荷载

转换为三角波的线性叠加后，运用基本的结构动力学知识即可解出体系位移情况。

本文重点内容：

分段函数的傅里叶展开
计算阻尼对振动的影响

二、理论推导

一般谐振荷载下的振动方程：

对周期荷载进行傅里叶展开：

其中的傅里叶系数：

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('傅里叶展开的外力函数图像');

结果