[转帖]静止的单自由度Duhamel积分
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%本程序只适用于初始状态为静止的单自由度Duhamel积分，内部采用Simpson积分
%%%%%%%可计算无阻尼和阻尼强迫震动
%%%%%%%输入可为数组或函数荷载
%%%%%%%只能得出位移时程图
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all;
%w=30
w=input('输入自振频率 ω: ');
%n=10
n=input('输入荷载步数n ：');
%m=96.6/32.3
m=input('输入单元质量m ：');
%T=0.05
T=input('输入外荷载持续时间 T：');
%si=0.0
si=input('输入阻尼系数ξ：');
%k=2700
k=input('输入刚度系数k：');
deta=T/n;
wD=w*(1-si^2)^0.5;
G=deta/(m*wD)/3;
t1=;
reply = input('输入荷载数组直接回车或敲击Y，输入函数荷载输入N： : ','s');
if isempty(reply)
reply = 'Y';
end
if reply=='Y'
p2=input('输入荷载数组并用[]抱住 为n＋1个元素：例如:\n');
%p2=
exp1=exp(-2*si*w*deta);
exp2=4*exp(-si*w*deta);
A(1)=p2(1)*cos(wD*0);
for i=2:n/2+1
A(i)=(A(i-1)+p2(2*i-3)*cos(wD*(2*deta*(i-2))))*exp1+p2(2*i-2)*cos(wD*(2*i*deta-3*deta))*exp2+p2(2*i-1)*cos(wD*2*(i-1)*deta);
end
B(1)=p2(1)*sin(wD*0);
for i=2:n/2+1
B(i)=(B(i-1)+p2(2*i-3)*sin(wD*(2*deta*(i-2))))*exp1+p2(2*i-2)*sin(wD*(2*i*deta-3*deta))*exp2+p2(2*i-1)*sin(wD*2*(i-1)*deta);
end
t1=;
disp('位移随时间的变化');
V=G*(A.*sin(wD*t1)-B.*cos(wD*t1))
disp('弹性力随时间的变化');
f=k*V
for i=1:n/2+1
p22(i)=p2(2*i-1);
end
disp('加速度随时间的变化');
v11=(p22-k*V)/m
subplot(1,2,1),plot(t1,V),title('位移'),grid on,
t2=;
sw=sin(w*t2);
cw=cos(w*t2);
disp('扩大四倍后位移随时间的变化');
y=G*A(n/2+1)*sw-G*B(n/2+1)*cw
subplot(1,2,2),
plot(t1,V,'r') ,hold on
plot(t2,y,'r');
title('加载时间增加四倍时的位移位移时程曲线')
grid on;
else
syms tao t;
y=input('函数荷载(关于未知数tao)：');
%y=-10000/8*tao+100;
v1=1/m/wD*int(y*exp(-si*w*(t-tao))*sin(wD*(t-tao)),tao,0,t);
v01=diff(v1,t);
v02=diff(v01,t)
t1=;
for i=1:n+1
v2(i)=subs(v1,t,t1(i));
v011(i)=subs(v01,t,t1(i));
v022(i)=subs(v02,t,t1(i));
end
disp('荷载作用时间内位移变化');
v2
disp('荷载作用时间内位移最大值');
max(v2)
disp('荷载作用时间内速度变化');
v011
disp('荷载作用时间内速度最大值');
max(v011)
disp('荷载作用时间内加速度变化');
v022
disp('荷载作用时间内加速度最大值');
max(v022)
subplot(2,3,1),plot(t1,v2),grid on,title('荷载作用时间内位移时程曲线');
subplot(2,3,2),plot(t1,v011),grid on,title('荷载作用时间内速度时程曲线');
subplot(2,3,3),plot(t1,v022),grid on,title('荷载作用时间内加速度时程曲线')
A=1/m/wD*int(y*exp(si*w*tao-si*w*t)*cos(wD*tao),tao,0,T);
B=1/m/wD*int(y*exp(si*w*tao-si*w*t)*sin(wD*tao),tao,0,T);
A1=subs(A,T);
B1=subs(B,T);
V33=A1*sin(wD*t)-B1*cos(wD*t);
V033=diff(V33,t);
V0033=diff(V033,t);
t2=;
v33=A1*sin(wD*t2)-B1*cos(wD*t2);
for i=1:4*n+1
V03(i)=subs(V033,t,t2(i));
V003(i)=subs(V0033,t,t2(i));
end
disp('扩大四倍后位移随时间的变化');
v33
disp('扩大四倍后位移最大值');
max(v33)
disp('扩大四倍后速度随时间的变化');
V03
disp('扩大四倍后速度最大值');
max(V03)
disp('扩大四倍后加速度随时间的变化');
V003
disp('扩大四倍后加速度最大值');
max(V003)
subplot(2,3,4),plot(t1,v2),hold on,title('加载时间增加四倍时的位移时程曲线'),plot(t2,v33),grid on;
subplot(2,3,5),plot(t1,v011),hold on,title('加载时间增加四倍时的速度时程曲线'),plot(t2,V03),grid on;
subplot(2,3,6),plot(t1,v022),hold on,title('加载时间增加四倍时的加速度时程曲线'),plot(t2,V003),grid on;
end
来自钢结构
[ 本帖最后由 suffer 于 2006-10-9 19:39 编辑 ] 谢谢，很好用，只是我想要个多自由度的。
页:
[1]