附件给出的是齿轮动力学振动模型和微分方程
其中x为要求的位移响应,Zeta为齿轮副相对阻尼系数,Zeta =0.07,m为等效质量m=3;Fn为法向载荷Fn=100N。
k(t)为齿轮副的时变啮合刚度
对k(t)我做了关于时间t的傅里叶变换:
k=750527161.016+1.2466e+008*sin(2*pi*f0*t-0.86045)+8.3314e+007*sin(4*pi*f0*t-0.2292)+4.9789e+007*sin(6*pi*f0*t+0.41062)+1.7784e+007*sin(8*pi*f0*t+1.0182)+7.633e+006*sin(10*pi*f0*t-1.2305);
这个方程是变系数微分方程,其中就k(t)一个时变系数,k(t)该如何表示?是否可把k(t)作为t的函数直接代入方程?
下面是我编的M文件和求解程序,各位大神帮我看看处理的对不对?谢谢
M文件
function Yd= fun7( t,x,k,zeta,m,Fn)
% FUN7.M: 时变微分方程例子
Fn=100;
zeta=0.07;
m=3;
n1=3000;
z1=45;
Tz=0.00044444;
f0=1/Tz;
k=750527161.016+1.2466e+008*sin(2*pi*f0*t-0.86045)+8.3314e+007*sin(4*pi*f0*t-0.2292)+4.9789e+007*sin(6*pi*f0*t+0.41062)+1.7784e+007*sin(8*pi*f0*t+1.0182)+7.633e+006*sin(10*pi*f0*t-1.2305);
dx=zeros(2,1);
dx(1)=x(2);
Yd=[x(2);Fn/m-2*zeta*sqrt(k/m)*x(2)-k/m*x(1)];
求解程序:
clear
m=3;
Tz=0.00044444;
f0=1/Tz;
tspan=linspace(0,0.00044444,120);
Y0=[1.5040e-007;0];
[t,XX]=ode45('fun7',tspan,Y0);
X=XX(:,1);
Y=XX(:,2);
plot(t,X),xlabel('经历时间/(s)'),ylabel('位移x'),figure;
plot(t,Y),xlabel('经历时间/(s)'),ylabel('速度v'),figure;
plot(X,Y)