clear all
clc
ebs = 0.05;%单自由度体系阻尼比
w = 10*pi;%单自由度系统无阻尼自振频率
w1 = w*(1-ebs^2)^0.5;%单自由度体统有阻尼自振频率计算公式
y0 = 0;%单自由度系统初始位移
y01 = 1;%单自由度系统初始速度
F = 300;%简谐荷载的干扰力的最大值
theta = 2*pi;%简谐荷载的频率
m =10;%单自由度体系质量
syms t
y = exp(-ebs*w*t)*(y0*cos(w1*t)+(y01+ebs*w*y0)/w1*sin(w1*t))+ ...
exp(-ebs*w*t)*theta*F/(m*((w^2-theta^2)^2+4*ebs^2*w^2*theta^2))*(2*ebs*w*cos(w1*t)+(2*ebs^2*w^2-(w^2-theta^2)/w1*sin(w1*t)))+...
F/(m*((w^2-theta^2)^2+4*ebs^2*w^2*theta^2))*((w^2-theta^2)*sin(theta*t)-2*ebs*w*theta*cos(theta*t));%结构力学单自由度在简谐荷载作用的激励下强迫振动的公式
dy = diff(y); %求一阶导数
d2y = diff(dy); %求二阶导数
y = matlabFunction(y);% 把变量表达式转化为句柄格式
dy = matlabFunction(dy);% 把变量表达式转化为句柄格式
d2y = matlabFunction(d2y);% 把变量表达式转化为句柄格式
T=3;%计算从0到3秒
dt = 0.001; %时间间隔
ti = (0:dt:T)'; %定义时间取值区间，将横向等差数列矩阵转置为竖向等差数列矩阵
y1 = y(ti); %计算y，即位移的值
y1=[y0 y1']';
y2 = dy(ti); %计算dy，即速度的值
y2=[y01 y2']';
y3 = d2y(ti); %计算d2y，即加速度的值
% y的画图
figure(1)
plot(ti,y1(1:length(ti)))
xlabel('t/s'),ylabel('y')%X坐标为时间，Y坐标为位移
% dy，即速度的画图
figure(2)
plot(ti,y2(1:length(ti)))
xlabel('t/s'),ylabel('dy')%X坐标为时间，Y坐标为速度
% d2y，即加速度的画图
figure(3)
plot(ti,y3)
xlabel('t/s'),ylabel('d2y')%X坐标为时间，Y坐标为加速度
format short g
disp(['时间，y,dy,d2y结果如下：'])
data = [ti,y1(1:length(ti)),y2(1:length(ti)),y3]%定义输出数据第一二三四列分别为时间，位移，速度，加速度
xlswrite('结构力学解答',data) %保存excel中，文件名为结构力学解答
评论0