0 个投票

I am trying to wrtie Shaft_speed as a function of time , so the shaft speed should increase to from 1 to 4000 as time goes on.
I also need the initial condition values X,Y and Z of the last speed to be the new intial condition of the new shaft speed. so far i have:
clc,clear,close all
J=3.37e-5; %inertia of rotor (kg/m2)
K1=0.7243; %Linear stiffness coef (N.m/rad)
K3=1409.7; %Nonlinear stiffness coef (N.m/rad3)
f=0.005;
Zeta=0.05;%damping ratio
Cmech=2*Zeta*(sqrt(K1/J))*J; %mechanical damping coef
Rint=82;
Rload=82;
ECF=sqrt(1.9e-6/Rint); %electromagnetic coupling factor
L=160e-3;%Inductance (mH)
%initial conditions
X = 0; %initial displacement
Y = 0; %initial Velocity
Z = 0; %initial .....
%shaft speed in RPM
Shaft_speed=0:1000;
%shaft speed in rad/s
Rad_speed=Shaft_speed.*((2*pi)/60);
%freq of harmonic fluctuations
w=2.*Rad_speed;
w_Hz=w./2*pi; %convert to Hertz
nt=1./w_Hz;
Time=0:nt/100:100*nt;
Iter=10;
for r=1:Iter
[t,q]=ode45(@duf,Time,[X Y Z]);
X=q(end,1);
Y=q(end,2);
Z=q(end,3);
Max_values(r) = max(q(:,1))
plot(t,q(:,1),'k')
xlabel('TIme,s')
ylabel('Relative Displacement, Rad')
grid
end
max=max(Max_values)
function qdot=duf(t,q) % q_ddot+(Cmech/J)*q_dot+(K1/J)*q+(K3/J)*x^3=f*(w^2)*cos(w*t)
J=3.37e-5; %inertia of rotor (kg/m2)
K1=0.7243; %Linear stiffness coef (N.m/rad)
K3=1409.7; %Nonlinear stiffness coef (N.m/rad3)
f=0.005;
Zeta=0.05;%damping ratio
Cmech=2*Zeta*(sqrt(K1/J))*J; %mechanical damping coef
Rint=82;
Rload=82;
ECF=sqrt(1.9e-6/Rint); %electromagnetic coupling factor
L=160e-3;%Inductance (mH)
%shaft speed in RPM
Shaft_speed=0:1000;
%shaft speed in rad/s
Rad_speed=Shaft_speed.*((2*pi)/60);
%freq of harmonic fluctuations
w=2.*Rad_speed;
w_Hz=w./2*pi; %convert to Hertz
qdot(1)=q(2);
qdot(2)=f.*w.^2.*cos(w.*t)-(Cmech/J).*q(2)-(K1/J).*q(1)-(K3/J).*q(1).^3- ECF.*q(3);
qdot(3)=(-ECF.*q(2)-(Rint+Rload).*q(3))./L;
qdot = [qdot(1);qdot(2);qdot(3)];
end
but i geet keep getting '-Error using odearguments (line 21)
When the first argument to ode45 is a function handle, the tspan argument must have at least two elements.
Error in ode45 (line 115)
odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);
Error in mainduffing (line 33)
[t,q]=ode45(@duf,Time,[X Y Z]);'
Thanks in advance

请先登录，再回答此问题。

 采纳的回答

VBBV
VBBV 2020-11-26

0 个投票

Shaft_speed=1; % specify a value instead of vector
%shaft speed in rad/s
Rad_speed=Shaft_speed*((2*pi)/60);
%freq of harmonic fluctuations
w=2*Rad_speed;
w_Hz=w/2*pi; %convert to Hertz
nt=1/w_Hz;
Time=0:nt:100*nt; % use coarse step length to reduce time of run
inside the for loop specify time as initial and final states
[t,q]=ode45(@duf,[0 max(Time)],[X Y Z]);
omit the element wise operation inside the function duf
qdot(2)=f*w^2*cos(w*t)-(Cmech/J)*q(2)-(K1/J).*q(1)-(K3/J)*q(1)^3- ECF*q(3);

6 个评论
显示 4更早的评论 隐藏 4更早的评论

VBBV
VBBV 2020-11-26
Inside the function duf you are trying to assign, different size dimensions for vectors, Since size of q is only 3 and shaft speed vector of different size,
Adedayo Odeleye
Adedayo Odeleye 2020-11-26
Thank you for the advice,
is there a way i can increase the value for Shaft_speed within the loop. So when the loop is finished calculating q at a certain speed, the final values of q are then the new initial condition for the next value of Shaft_speed
VBBV
VBBV 2020-11-26
Shaft_speed=0.2:1000; % specify a arbitrary small value for shaft speed, bcoz a 0 speed will produce infinite time value
for i = 1:length(Shaft_speed)
%shaft speed in rad/s
S_Speed = Shaft_speed(i);
Rad_speed=Shaft_speed(i)*((2*pi)/60);
%freq of harmonic fluctuations
w=2*Rad_speed;
w_Hz=w/2*pi; %convert to Hertz
nt=1/w_Hz;
Time=0:nt:100*nt;
Iter=10;
for r=1:Iter
[t,q]=ode45(@duf,[0 max(Time)],[X Y Z]);
X=q(1); % Final values are assigned to initial values for next iteration
Y=q(2);
Z=q(3);
Max_values(r) = max(q(:,1))
plot(t,q(:,1),'k')
xlabel('TIme,s')
ylabel('Relative Displacement, Rad')
grid
end
end
function qdot=duf(t,q,S_Speed) % q_ddot+(Cmech/J)*q_dot+(K1/J)*q+(K3/J)*x^3=f*(w^2)*cos(w*t)
...
...
end
Pass the shaft speed as function argument
Adedayo Odeleye
Adedayo Odeleye 2020-11-26
Thank you very much! Really apericiate your help. The script is running fine.
Just one last question, what is the implication of having a different value for Shaft_speed in the function because am not able to mirror 'Shaft_speed=0.2:1000' into the function part, I can only write a single value unless I get error:
'Error using ^
One argument must be a square matrix and the other must be a scalar. Use POWER (.^) for elementwise power.
Error in Model>duf (line 65)
qdot(2)=f*w^2*cos(w*t)-(Cmech/J)*q(2)-(K1/J)*q(1)-(K3/J)*q(1)^3- ECF/J*q(3);
Error in odearguments (line 90)
f0 = feval(ode,t0,y0,args{:}); % ODE15I sets args{1} to yp0.
Error in ode45 (line 115)
odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);
Error in Model (line 32)
[t,q]=ode45(@duf,[0 max(Time)],[X Y Z]);'
VBBV
VBBV 2020-11-26
Since you are using Shaft_speed as vector inside duf function comment that line. And use S_Speed in place of Shaft_speed.
If you use it as vector inside the function the matrix dimensions do not match with output vector q size I.e.3
VBBV
VBBV 2020-11-26
Instead pass shaft spped as function argument I.e. S_Speed
Also define a coarse step size for Shaft_speed in the beginning to reduce computation time e.g.
%if true
Shaft_speed = 0.2:200:1000;

请先登录，再进行评论。

更多回答（1 个）

Walter Roberson
Walter Roberson 2020-11-26

0 个投票

Shaft_speed=0:1000;
That is a vector, with first element 0.
%shaft speed in rad/s
Rad_speed=Shaft_speed.*((2*pi)/60);
So that is a vector, with first element 0.
%freq of harmonic fluctuations
w=2.*Rad_speed;
So w is a vector, with first element 0.
w_Hz=w./2*pi; %convert to Hertz
so w_Hz is a vector, with first element 0.
nt=1./w_Hz;
so nt is a vector, with first element 1/0 -> inf
Time=0:nt/100:100*nt;
and the right hand side is then 0:nt(1):100*nt(1) -> 0:inf/100:100*inf -> 0:inf:inf which is just the element 0 .
When you use a vector in a colon operation, the value used is the first one, which happens to be inf in this case.

1 个评论
显示 -1更早的评论 隐藏 -1更早的评论

Adedayo Odeleye
Adedayo Odeleye 2020-11-26
Thank you for the help, i have replaced Shaft_speed with 1:1000 now i am getting :
In an assignment A(:) = B, the number of elements in A and B must be the same.
Error in mainduffing>duf (line 72)
qdot(2)=f.*w.^2.*cos(w.*t)-(Cmech/J).*q(2)-(K1/J).*q(1)-(K3/J).*q(1).^3- (ECF/J).*q(3);
Error in odearguments (line 90)
f0 = feval(ode,t0,y0,args{:}); % ODE15I sets args{1} to yp0.
Error in ode45 (line 115)
odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);
Error in mainduffing (line 32)
[t,q]=ode45(@duf,[Time],[X Y Z]);

请先登录，再进行评论。

类别

在 帮助中心 和 File Exchange 中查找有关 Work with Components 的更多信息

产品

版本

R2017a

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by