%状态方程
%给程序中的参数赋值
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
l=1;b=0.05;h=0.05;%梁的尺寸
a1=523.29;a2=1.868*10^7;a3=2.186*10^9;%本构方程中的材料参数
tm=15;%相变温度
i1=b*h^3/12;i2=b*h^5/80;i3=b*h^7/448;%梁的截面矩
%给无量纲量赋值
arf2=a2/(a1*tm);arf3=a3/(a1*tm);
c2=arf2*i2/(i1*l^2);c3=arf3*i3/(i1*l^4);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%ctd无量纲温度,v无量纲速度
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
global ctd v k d q;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%b1,b2,b3为状态方程中的系数
%%%%%%%%%%%%%%%%%%%%%%%%%
b1=pi^4*(ctd-1)-pi^2*v^2;
b2=3/4*pi^8*c2;
b3=5/8*pi^12*c3;
%%%%%%%%%%%%%%%%%%%%%%%%%
dy=zeros(3,1);
dy(1)=y(2);
dy(2)=-k*y(2)-b1*y(1)+b2*y(1)^3-b3*y(1)^5+d*sin(y(3));
dy(3)=q;
%相图及时域图
%ctd无量纲温度,v无量纲速度,k无量纲阻尼比,d无量纲激振力幅值,q无量纲激振力频率
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
global ctd v k d q;
ctd=1.7;v=0;k=0.05;d=0.02;q=0.5;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
tspan=[0:0.01*2*pi:5000];
y0=[0.001,0,0];
[t,y]=ode45('movingbeamforce',tspan,y0);
figure(1)
hold on
plot(y(end-1000:end,1),y(end-1000:end,2))
xlabel('无量纲位移'),ylabel('无量纲速度')
hold on
figure(2)
hold on
plot(t(end-1000:end),y(end-1000:end,1))
xlabel('无量纲时间'),ylabel('无量纲位移')
hold on
figure(3)
hold on
plot(t(end-1000:end),y(end-1000:end,2))
xlabel('无量纲时间'),ylabel('无量纲速度')
hold on
%Poincare截面程序
function mbeamforcepoincare
global ctd v k d q;
ctd=1.7;v=0;k=0.05;d=0.02;q=0.5;
y0=[0.001,0,0];
tspan=[0:0.01*2*pi:5000];
[t,y]=ode45('movingbeamforce',tspan,y0);
m=zeros(790,1);
n=zeros(790,1);
for i=1:790
m(i,1)=y(100*i,1);
n(i,1)=y(100*i,2);
end
plot(m,n,'k.','markersize',1); |