How to fix error in my DAE program?

2 次查看（过去 30 天）
Siddharth Jain
Siddharth Jain 2023-1-11
评论： Torsten 2023-1-11
The follwing is my main code- I am getting error -
Not enough input arguments.
Error in @(t,y,ydot)reduced21(t,y,ydot,T_a)
Error in odearguments (line 90)
f0 = feval(ode,t0,y0,args{:}); % ODE15I sets args{1} to yp0.
Error in ode15s (line 152)
odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);
function [yp] = reduced21(t,y,ydot,T_a)
% Define the initial state vector y0 and the initial derivative vector ydot0
y0 = zeros(12,1);
ydot0 = zeros(12,1);
beta = 13*(pi/180); %Helix angle (rad)
speed = 1000/60;
Freq = 1000*20/60;
% Common parameters
theta_a_vec = -speed*2*pi*t; %torsional angle of driver gear (rad)
% theta_p = 22.9; %torsional angle of driven gear (rad)
theta_p = pi/6 + speed*2*pi*t;
%e = 20e-6; %circumferential relative displacement of the teeth (m)
e = 0;
z = e*tan(beta); %axial tooth displacement caused by internal excitation (m)
e_a = 48e-6; %circumferential displacement of the driver gear (m)
e_p = 48e-6; %circumferential displacement of the driver gear (m)
ks = 0.9e8 + 20*sin(2*pi*Freq*t); %Shear stiffness
Cs = 0.1; %Shear damping
% Driver gear
m_a = 0.5; %mass of driver gear (kg)
c_ay=500; %Damping of driver gear in y direction (Ns/m)
c_az=500; %Damping of driver gear in z direction (Ns/m)
k_ay= 8e7; %Stiffness of driver gear in y direction (N/m)
k_az= 5e7; %Stiffness of driver gear in z direction (N/m)
R_a = 51.19e-3; %Radius
I_a = 0.5*m_a*(R_a^2); %Moment of inertia of driver gear (Kg.m3)
y_ac= e_a + theta_a_vec*R_a; %circumferential displacement at the meshing point on the driver gear (m)
y_pc = e_p - theta_a_vec*R_a; %circumferential displacement at the meshing point on the driven gear (m)
z_a = e_a*tan(beta);
z_p = e_p*tan(beta);
z_ac = (z_a-y_ac)*tan(beta); %axial displacements of the meshing point on the driving gear (m)
z_pc = (z_p+y_pc)*tan(beta); %axial displacements of the meshing point on the driven gear (m)
% Excitation forces
Fy=ks*cos(beta)*(y_ac-y_pc-e) + Cs*cos(beta)*2*R_a*speed*2*pi; %axial dynamic excitation force at the meshing point (N)
Fz=ks*sin(beta)*(z_ac-z_pc-z) - Cs*sin(beta)*2*tan(beta)*R_a*speed*2*pi; %circumferential dynamic excitation force at the meshing point (N)
%Time interpolation - Needed for solution of equations (It basically uses
%your torque and prescribed time matrices to generate a time matrix to be
%used in the ODE solver)
time = 0:0.00001:0.6;
Torque = interp1(time,T_a,t);
Opp_Torque = 68.853; % Average torque value
%Driver gear equations
yp = zeros(12,1); %vector of 12 equations
yp(1) = y(2);
yp(2) = (-Fy - k_ay*y(1) - c_ay*y(2))/m_a; %acceleration in y (up and down) direction (m/s2)
yp(3) = y(4);
yp(4) = (-Fz - k_az*y(3) - c_az*y(4))/m_a; %acceleration in z (side to side) direction (m/s2)
yp(5) = y(6);
yp(6) = (-Torque - Fy*R_a + Opp_Torque)/I_a; % angular acceleration (rad/s2)
% Driven gear
m_p = 0.5; %mass of driven gear (kg)
c_py=500; %Damping of driven gear in y direction (Ns/m)
c_pz=500; %Damping of driven gear in z direction (Ns/m)
k_py= 9.5e7; %Stiffness of driven gear in y direction (N/m)
k_pz =9e7; %Stiffness of driven gear in z direction (N/m)
R_p = 139.763e-3; %Radius
I_p = 0.5*m_p*(R_p^2); %Moment of inertia of driven gear (Kg.m3)
n_a = 20; % no of driver gear teeth
n_p = 60; % no of driven gear teeth
i = n_p/n_a; % gear ratio
%Driven gear equations
yp(7) = y(8);
yp(8) = (-Fy - k_py*y(7) - c_py*y(8))/m_p; %acceleration in y (up and down) direction (m/s2)
yp(9) = y(10);
yp(10) = (-Fz - k_pz*y(9) - c_pz*y(10))/m_p; %acceleration in y (side to side) direction (m/s2)
yp(11) = y(12);
yp(12) = (-Torque*i - Fy*R_p + Opp_Torque*i)/I_p; % angular accelration (rad/s2)
y0 = [0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0]; % matrix of intial values of y
ydot0 = [0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0]; %matrix of initial vlaues of ydot0
end
My command window code-
tic
A = load("torque_for_Sid_60s.mat");
T_a = A.torque_60s(1:60001); %Torque (Nm)
speed = 1000/60;
tspan=0:0.00001:0.6;
%tspan = [0 12];
y0 = zeros(12,1);
ydot0 = zeros(12,1);
y0 = zeros(12,1);
[t,y,ydot] = ode15s(@(t,y,ydot) reduced21(t,y,ydot,T_a),tspan,y0,ydot0);
toc
% Driver gear graphs
figure
subplot(3,1,1)
plot(t,y(:,2))
ylabel('(y) up and down acceleration (m/s2)')
xlabel('Time')
title('acceleration in y direction vs time')
grid on
subplot(3,1,2)
plot(t,y(:,4))
ylabel('(z) side to side acceleration (m/s2)')
xlabel('Time')
title('acceleration in z direction vs time')
grid on
subplot(3,1,3)
plot(t,y(:,6))
ylabel('angular acceleration (rad/s2)')
xlabel('Time')
title('angular acceleration vs time')
grid on
% Driven gear graphs
figure
subplot(3,1,1)
plot(t,y(:,8))
ylabel('(y) up and down acceleration (m/s2)')
xlabel('Time')
title('acceleration in y direction vs time')
grid on
subplot(3,1,2)
plot(t,y(:,10))
ylab
el('(z) side to side acceleration (m/s2)')
xlabel('Time')
title('acceleration in z direction vs time')
grid on
subplot(3,1,3)
plot(t,y(:,12))
ylabel('angular acceleration (rad/s2)')
xlabel('Time')
title('angular acceleration vs time')
grid on
% Torque graph
figure
plot(t,T_a(1:60001))
ylabel('Torque (Nm)')
xlabel('Time (sec)')
title('Torque vs time')
grid on

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

回答（1 个）

Torsten
Torsten 2023-1-11
编辑：Torsten 2023-1-11
The call to ode15s is incorrect:
[t,y] = ode15s(@(t,y) reduced21(t,y,T_a),tspan,y0);
instead of
[t,y,ydot] = ode15s(@(t,y,ydot) reduced21(t,y,ydot,T_a),tspan,y0,ydot0);
and the function header has to be modified accordingly:
function yp = reduced21(t,y,T_a)
...
end
  4 个评论
显示 2更早的评论
Siddharth Jain
Siddharth Jain 2023-1-11
ydot0 is a vector of differential of the initial values in y0
Torsten
Torsten 2023-1-11
It's not needed (and you cannot prescribe it).
ydot0 will be the vector yp that is calculated in reduced21 by inserting (t,y) = (0,y0) into the equations.

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Programming 的更多信息

标签

产品

Community Treasure Hunt

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

Start Hunting!

Translated by