核心是解一个三自由度微分方程组,解出θ关于时间t的变化图像.设定一个变转速的工况,图中g(t)是功率随时间t变化函数,wm是根据功率求得的角速度.图一是主程序,图二是参数和方程,图三是求解函数,系统提示我无法将表达式转换为双数组,想问一下我这个错在哪里有没有修改办法或者建议,万分感谢
clear;clc;close all;
%初始参数，该参数变化，myfun程序要对应调整
syms t;
g(t)=(1460*9.8+1.04*1460*2.67+0.32*2.25*(10+2.67*t)^2/21.15)*(10+2.67*t)/3600;
wm=(g(t)*(10+2.67*t)*9550/(3600*0.4064/(0.93*14)))*pi*2/60;
zi=4;%内转子极对数
zo=14;%外转子极对数
z1=32;%小齿轮齿数
z2=99;%大齿轮齿数
ms=4.62e-3;
% r2=ms*z1/2;%小齿轮分度圆
% r3=ms*z2/2;%大齿轮分度圆
I=z2/z1;%机械齿轮速比
% G=zo/zi;
% ms1=4.12e-3;
G=(zo+zi)/zi;%磁齿轮速比
%一、计算过程
tmax=10;
% tmax=900*2*pi/wm;% 计算时长，按输入轴的周期数计算
bc=2e-5;
y0=[0;0;0;wm;wm/G;wm/G/I];
% aa=0:1e-4:tmax;
% [t,y]=ode45(@myfun3, aa, y0);
[t,y]=odeRK4sys('myfun3',tmax,bc,y0);
figure;plot(t,y(:,1)/G-y(:,2));
title('时域响应图');
xlabel('时间 (s)');
ylabel('扭转振幅(rad)');%磁齿轮的相对弹性角度
figure;plot(t,y(:,2)/I-y(:,3));
title('时域响应图');
xlabel('时间 (s)');
ylabel('扭转振幅(rad)');%机械齿轮的相对弹性转角;
figure;plot(t,y(:,4),t,y(:,5),t,y(:,6));
title('时域响应图');
xlabel('时间 (s)');
ylabel('转速(rad/s)');
function dydt=myfun3(t,y)
syms t;
g(t)=(1460*9.8+1.04*1460*2.67+0.32*2.25*(10+2.67*t)^2/21.15)*(10+2.67*t)/3600;
wm=(g(t)*(10+2.67*t)*9550/(3600*0.4064/(0.93*14)))*pi*2/60;
zi=4;%内转子极对数
zo=14;%外转子极对数
z1=32;%小齿轮齿数
z2=99;%大齿轮齿数
arf0=20/180*pi;%压力角
ms=4.62e-3;
r2=ms*z1/2;%小齿轮分度圆
r3=ms*z2/2;%大齿轮分度圆
rb2=r2*cos(arf0);%小齿轮基圆
rb3=r3*cos(arf0);%大齿轮基圆
I=z2/z1;%机械齿轮速比
G=(zo+zi)/zi;%磁齿轮速比
Tm=(1460*9.8+1.04*1460*2.67+0.32*2.25*(10+2.67*t)^2/21.15)*0.4064/(0.93*14);%电机输入扭矩
TL=Tm*I*G;%变速器负载
%3 转动惯量
IM=0.081197;%电机转子转动惯量kgm2
I0=0.181197;%磁齿轮低速级转动惯量kgm2
I1=0.013663;%磁齿轮高速级转动惯量kgm2
I2=0.0230596;%机械小齿轮转动惯量
I3=1.9825692;%机械大齿轮转动惯量
IL=116;%整车等效转动惯量
%4 刚度与阻尼
%齿轮的阻尼和刚度
w1=wm*G;%机械小齿轮输入转速
kp=10.51e8;
kb=0.25e8;%两个齿轮的刚度曲线为正弦曲线，
km=kp+kb*sin(w1*t*z1);
sg=0.16;%齿轮阻尼比
cm=2*sg*sqrt(kp*I2*I3/(I2*rb3^2+I3*rb2^2));
b=1e-4;
xx=rb2*y(2)-rb3*y(3);
dxx=rb2*y(5)-rb3*y(6);
if xx-b>0
fx=xx-b;
elseif xx+b<0
fx=xx+b;
else
fx=0;
end
F23=km*fx+cm*dxx;
T0=1.7*Tm*G;
Tc1=0.1*T0*sin(wm*t);
Tc=0.1*(zi*y(4)-(zi+zo)*y(5));
T21=T0*sin(zi*y(1)-(zi+zo)*y(2))+Tc+Tc1;
y(7)=(Tm-T21/G)/(I0+IM);
y(8)=(T21-r2*F23)/(I2+I1);
y(9)=(-TL+r3*F23)/(I3+IL);
dydt=[y(4)
y(5)
y(6)
y(7)
y(8)
y(9)];
end
function [t,y] = odeRK4sys(diffeq,tn,h,y0)
t = (0:h:tn)'; % Column vector of elements with spacing h
nt = length(t); % number of steps
neq = length(y0); % number of equations that are simultaneously advanced
y = zeros(nt,neq); % Preallocate y for speed
y(1,:) = y0(:)'; % Assign IC. y0(:) is column, y0(:)' is row vector
h2 = h/2; h3 = h/3; h6 = h/6; % Avoid repeated evaluation of constants
k1 = zeros(neq,1); k2 = k1; % Preallocate memory for the Runge-Kutta
k3 = k1; k4 = k1; ytemp = k1; % coefficients and a temporary vector
% Outer loop for all steps: j = time step index; n = equation number index
for j=2:nt
told = t(j-1); yold = y(j-1,:)'; % Temp variables, yold is column vector
k1 = feval(diffeq,told,yold); % Slopes at the starting point
for n=1:neq
ytemp(n) = yold(n) + h2*k1(n); % Estimate all y's at midpoint
end
k2 = feval(diffeq,told+h2,ytemp); % 1st estimate of slopes at midpoint
for n=1:neq
ytemp(n) = yold(n) + h2*k2(n); % 2nd estimate of all y's at midpoint
end
k3 = feval(diffeq,told+h2,ytemp); % 2nd estimate of slopes at midpoint
for n=1:neq
ytemp(n) = yold(n) + h*k3(n); % Estimate y at end point
end
k4 = feval(diffeq,told+h,ytemp); % Estimate of slopes at endpoint
for n=1:neq % Simultaneously advance all y's
y(j,n) = yold(n) + h6*(k1(n)+k4(n)) + h3*(k2(n)+k3(n));
end
end
