How to convert ode45 to Differential Algebraic Equations?
2 次查看(过去 30 天)
显示 更早的评论
I have 12 equations with 6 DOF ( see yp) that i want to convert to Differential Algebraic Equation format. How can I do this?
function [yp] = reduced(t,y,T_a)
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)
e = 20e-6; %circumferential relative displacement of the teeth (m)
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 = 5.3*e-2; %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: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 = 2.5*e-1; %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= 8e7; %Stiffness of driven gear in y direction (N/m)
k_pz =5e7; %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
Opp_Torque = 68.853; % Average torque value
%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)
end
3 个评论
Torsten
2022-12-14
Solving differential equations and algebaic equations simultaneously is not possible with ODE45. If this is what your question aims at, you will have to use ODE15S or ODE23T, as John D'Errico colourfully explained.
回答(1 个)
John D'Errico
2022-12-14
编辑:John D'Errico
2022-12-14
I would like to fly to Los Angeles. However, I don't want to take a plane, and since I have often taken a train to destinations, does anyone know how to make a train fly?
The point is, a train is not designed to fly. No matter how hard you try, you will not get it off the tracks without rather unhappy consequences.
In your case, you have a DAE: a Differential Algebraic Equation system. So you use a tool designed to solve a DAE. As I recall, ODE15S will do that, and ODE23T. Maybe some others, but those are the ones that spring to mind.
help ode15s
help ode23t
If you want to fly, you need to use something designed to fly.
0 个评论
另请参阅
产品
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!