使用ODE45求解齿轮系统动力学方程后结果发散

大小的小

1.1参数定义及动力学方程降阶
# P* B$ M! X) }% z$ P# h9 O" ffunction [dx,ff1,ff2]=myfun(t,x)
* z. k$ t2 A3 [8 }; Rt
beita=26;
mn=0.004;
z1=46;                       
z2=43;                       
z3=122;                     
T_in=200;
T_out=80;
3 u3 P* d' g' S1 Droug1=7.8E3;            
% m4 A8 a6 r' j& r  Xroug2=7.8E3;            
5 Y5 j( N* l# A7 s$ Groug3=7.8E3;            
. t8 P9 u* Z  u( I4 U0 l, Lalphan=20;                                                   
1 g6 {; b2 o- T( S+ }alphat=atand(tand(alphan)/cosd(beita));      
d1=z1*mn/cosd(beita)/1000;                  
db1=d1*cosd(alphat)/1000;                    
5 o- y' L" Y# {1 q5 T2 c! i6 V7 dd2=z2*mn/cosd(beita)/1000;                  
8 P, N* U/ z7 g1 Y# y- Y, V& }- Mdb2=d2*cosd(alphat)/1000;                    
d3=z3*mn/cosd(beita)/1000;                  
db3=d3*cosd(alphat)/1000;                                 
+ S# f" v- a2 U& F6 pbp1=116/1000;                                       
7 {  R8 I% @1 _) Dbp2=116/1000;                                       
bp3=116/1000;                                       
bp=116/1000;
I1=((roug1*pi*(d1/2)^2*bp1)*(d1/2)^2)/2;  
4 q: u4 e* x: N; d4 u' _I2=((roug2*pi*(d2/2)^2*bp2)*(d2/2)^2)/2;  
9 n9 r  V9 A0 u/ n$ UI3=((roug3*pi*(d3/2)^2*bp3)*(d1/2)^2)/2;  
4 K4 j" Z( N/ ]- \: Gm1=roug1*pi*(d1/2)^2*bp1;                                
m2=roug2*pi*(d2/2)^2*bp2;                              
# }# O2 x+ j- }/ f" Qm3=roug3*pi*((d3)-(d1+d2))^2*bp3;                  
# M6 G: s; P; qr1=d1*cosd(alphat)/2;                                          
r2=d2*cosd(alphat)/2;                                          
6 d) a( f9 k& b7 |5 }. }; er3=d3*cosd(alphat)/2;                                          
$ W. d% `/ p+ J# N/ L; K5 ]  S( p7 ufai_sp1x=90;
fai_sp1y=0;
0 T: Z2 B! Y  Q9 q6 ifai_p1rx=-130;
- j3 n$ _8 [$ c' H6 ifai_p1ry=-220;
kesaiz=0.05;
8 {* t# X2 `; y/ bkesain=0.07;
kp1x=1e8;
kp1y=kp1x;
) \0 |0 J' x- Mcp1x=2*kesaiz*((kp1x*m2)^0.5);
8 u: y6 D1 A% z! y! ^cp1y=2*kesaiz*((kp1y*m2)^0.5);
% j' o! y3 X9 U0 P/ l/ tksx=1e8;
ksy=1e8;                                   
9 d: K0 t0 q1 Z6 Rcsx=2*kesaiz*((ksx*m1)^0.5);
csy=2*kesaiz*((ksy*m1)^0.5);
. h6 _% e+ `- n3 ?+ Ikrx=1e8;
kry=krx;
% h2 s3 {$ S, L2 Scrx=2*kesaiz*((krx*m3)^0.5);   
cry=crx;
Tmesh=2*pi/z1;
kp1r =1e6;
! V2 D  M, `) @& F( C) [csp1=2*0.07*((ksp1/(1/m1+1/m2))^0.5);
cp1r=2*0.07*((ksp1/(1/m2+1/m3))^0.5);                 
& \( o1 K* v$ M0 @$ Q- t% t%%%%%%%%%%%%%%%%%%%%%%%%%%%%以上为参数定义，可忽略不看，谢谢！
& r' u- k+ J  iesp1=1e-6;
ep1r=1e-6;
delta_sp1=((x(1)-x(7))*cosd(fai_sp1x)+(x(3)-x(9))*cosd(fai_sp1y)+r1*x(5)+r2*x(11))*cosd(beita)+esp1;
delta_sp11=((x(2)-x(8))*cosd(fai_sp1x)+(x(4)-x(10))*cosd(fai_sp1y)+r1*x(6)+r2*x(12))*cosd(beita)+esp1;
9 m! j7 w. {/ h0 r, Hdelta_p1r=((x(7)-x(13))*cosd(fai_p1rx)+(x(9)-x(15))*cosd(fai_p1ry)+r3*x(17)-r2*x(11))*sind(beita)-ep1r;
& ~+ H+ O' j1 |delta_p11r=((x(8)-x(14))*cosd(fai_p1rx)+(x(10)-x(16))*cosd(fai_p1ry)+r3*x(18)-r2*x(12))*sind(beita)-ep1r;
% \4 F( y- d$ _' W%%%%%%%%%%动力学方程
3 M: n9 m9 W2 Q. Wdx=zeros(18,1);
dx(1)=x(2);
dx(2)=(1500*cosd(fai_sp1x)-(csx*x(2))-ksx*x(1)-(csp1*delta_sp11+ksp1*delta_sp1)*cosd(beita)*cosd(fai_sp1x))/m1;
6 o3 F) ^4 `4 Tdx(3)=x(4);
dx(4)=(1500*cosd(fai_sp1y)-(csx*x(4))-ksy*x(2)-(csp1*delta_sp11+ksp1*delta_sp1)*cosd(beita)*cosd(fai_sp1y))/m1;
dx(5)=x(6);
" k8 i+ K! h* adx(6)=(400-(csp1*delta_sp11+ksp1*delta_sp1)*cosd(beita)*r1+T_in)/I1;
dx(7)=x(8);
# ~; S$ t9 {& Hdx(8)=((-300*cosd(fai_sp1x)+200*cosd(fai_p1rx)-(cp1x*x(8))-kp1x*x(7)+(csp1*delta_sp11+ksp1*delta_sp1)*cosd(beita)*cosd(fai_sp1x)-(cp1r*delta_p11r+kp1r*delta_p1r))*cosd(beita)*cosd(fai_p1rx))/m2;                %第一个行星齿轮x方向
dx(9)=x(10);
& L$ r, B  k9 ^dx(10)=((-300*cosd(fai_sp1y)+200*cosd(fai_p1ry)-(cp1y*x(10))-kp1y*x(9)+(csp1*delta_sp11+ksp1*delta_sp1)*cosd(beita)*cosd(fai_sp1y)+(cp1r*delta_p11r+kp1r*delta_p1r))*cosd(beita)*cosd(fai_p1ry))/m2;                 %第一个行星齿轮y方向
0 [8 r  t$ L9 {" s" N% Ydx(11)=x(12);
dx(12)=(120-((csp1*delta_sp11+ksp1*delta_sp1)+(cp1r*delta_p11r+kp1r*delta_p1r))*cosd(beita)*r2)/I2;
. F  f" {$ L" q# V, }' {dx(13)=x(14);
dx(14)=((-200*cosd(fai_p1rx)-(crx*x(14))-krx*x(13)+(cp1r*delta_p11r+kp1r*delta_p1r))*cosd(beita)*cosd(fai_p1rx))/m3;
dx(15)=x(16);
dx(16)=((-200*cosd(fai_p1ry)-(crx*x(16))-krx*x(15)+(cp1r*delta_p11r+kp1r*delta_p1r))*cosd(beita)*cosd(fai_p1ry))/m3;
# u! U4 G* x6 Zdx(17)=x(18);
2 \( q  D8 K8 @0 a) s, n0 W4 w" Zdx(18)=(80-(cp1r*delta_p11r+kp1r*delta_p1r)*cosd(beita)*r3-T_out)/I3;

  L( u0 |, `! j* \2 U' p

1.2 ode程序
clc;
clear all
; H% o2 K* T* _0 _% c7 b) Lx0=zeros(18 ,1)
3 _; X9 ^2 ?9 h/ H) ~, [( ]- x[t,x] = ode45('myfun',[0:0.0001:10],x0);
* [' ]9 B8 V! I  U. m# X/ efigure
plot(t,x(:,1))


3.绘图结果如下，为什么画出来是一条直线，而且图中结果没有计算到规定的时间



5 n( O( N3 c7 r9 [1 P7 k
# j0 L) \, w2 `$ M( h" E3 ]

Uifhjvv

帮你顶一下

IBB-EUT

你这个自定义函数不完整吧，找不到对 ksp1 这个变量的有效赋值，csp1等等需要调用ksp1的赋值也没法做

shuddkk

function [dx,ff1,ff2]=myfun(t,x)
1 \% n& O/ ^! t& K8 `9 D  lt
) |6 j. F: @0 d: H4 Z: Abeita=26;
mn=0.004;
4 E% F* b6 l" C4 B$ \1 Kz1=46;
( V4 ~) S4 M, h/ B8 az2=43;
9 N9 a6 K& Q) ez3=122;
T_in=200;
T_out=80;
roug1=7.8E3;
roug2=7.8E3;
roug3=7.8E3;
alphan=20;
alphat=atand(tand(alphan)/cosd(beita));
! L5 z( z9 M  p" b0 L7 id1=z1*mn/cosd(beita)/1000;
db1=d1*cosd(alphat)/1000;
d2=z2*mn/cosd(beita)/1000;
6 q* g( j9 ]& f3 c4 R2 P9 s0 @db2=d2*cosd(alphat)/1000;
d3=z3*mn/cosd(beita)/1000;
db3=d3*cosd(alphat)/1000;
4 L4 \. ?# O& ?bp1=116/1000;
bp2=116/1000;
, N0 c( G+ P8 R6 y* n5 Jbp3=116/1000;
3 c( K" ^" a( Q# y! p6 e- m- mbp=116/1000;
, s' k- P. g9 w& g) U4 H7 G4 @I1=((roug1*pi*(d1/2)^2*bp1)*(d1/2)^2)/2;
I2=((roug2*pi*(d2/2)^2*bp2)*(d2/2)^2)/2;
I3=((roug3*pi*(d3/2)^2*bp3)*(d1/2)^2)/2;
m1=roug1*pi*(d1/2)^2*bp1;
+ ?& W# O: ]( \m2=roug2*pi*(d2/2)^2*bp2;
m3=roug3*pi*((d3)-(d1+d2))^2*bp3;
- b) k. K3 e( _% n' e  S& Nr1=d1*cosd(alphat)/2;
r2=d2*cosd(alphat)/2;
$ W( S0 z7 t. C% t* R2 e, gr3=d3*cosd(alphat)/2;
$ H+ r) P3 x& c) Kfai_sp1x=90;
' H( w9 X0 ]" d' F" C: C, ufai_sp1y=0;
9 {2 o% t6 P' Ffai_p1rx=-130;
fai_p1ry=-220;
kesaiz=0.05;
kesain=0.07;
ksp1 =1e6;
) }* J8 v: N: s. H9 R1 fkp1r =1e6;
kp1x=1e8;
kp1y=kp1x;
7 m  H, `+ _. r* F# E% k3 vcp1x=2*kesaiz*((kp1x*m2)^0.5);
cp1y=2*kesaiz*((kp1y*m2)^0.5);
ksx=1e8;
5 \& U, V; L6 r& `9 zksy=1e8;
' t3 e: U* o4 |6 Hcsx=2*kesaiz*((ksx*m1)^0.5);
csy=2*kesaiz*((ksy*m1)^0.5);
krx=1e8;
) D# S& K& u+ [5 ^- y& Z2 lkry=krx;
2 J* y% I' |4 U4 X0 ~. U% Ocrx=2*kesaiz*((krx*m3)^0.5);
8 J( N3 z/ y3 q0 k) z5 xcry=crx;
Tmesh=2*pi/z1;
$ S9 E, i% y4 B3 p0 ykp1r =1e6;
, W4 l: O  _* m( Lcsp1=2*0.07*((ksp1/(1/m1+1/m2))^0.5);
cp1r=2*0.07*((ksp1/(1/m2+1/m3))^0.5);
7 h9 ~$ [4 e3 V# ]1 r; Q  M%%%%%%%%%%%%%%%%%%%%%%%%%%%%以上为参数定义，可忽略不看，谢谢！
esp1=1e-6;
. k. g4 w/ K* z0 v. E  zep1r=1e-6;
delta_sp1=((x(1)-x(7))*cosd(fai_sp1x)+(x(3)-x(9))*cosd(fai_sp1y)+r1*x(5)+r2*x(11))*cosd(beita)+esp1;
5 B7 D( K0 A+ ^delta_sp11=((x(2)-x(8))*cosd(fai_sp1x)+(x(4)-x(10))*cosd(fai_sp1y)+r1*x(6)+r2*x(12))*cosd(beita)+esp1;
% Y  s- I, n  M0 Y" b: c, idelta_p1r=((x(7)-x(13))*cosd(fai_p1rx)+(x(9)-x(15))*cosd(fai_p1ry)+r3*x(17)-r2*x(11))*sind(beita)-ep1r;
delta_p11r=((x(8)-x(14))*cosd(fai_p1rx)+(x(10)-x(16))*cosd(fai_p1ry)+r3*x(18)-r2*x(12))*sind(beita)-ep1r;
%%%%%%%%%%动力学方程
dx=zeros(18,1);
5 V: R' S  d; {dx(1)=x(2);
) {" Y: a+ E1 U- `. n- gdx(2)=(1500*cosd(fai_sp1x)-(csx*x(2))-ksx*x(1)-(csp1*delta_sp11+ksp1*delta_sp1)*cosd(beita)*cosd(fai_sp1x))/m1;
dx(3)=x(4);
$ l9 H3 A+ D; ?, U) E2 X" ydx(4)=(1500*cosd(fai_sp1y)-(csx*x(4))-ksy*x(2)-(csp1*delta_sp11+ksp1*delta_sp1)*cosd(beita)*cosd(fai_sp1y))/m1;
+ [% v) @7 c) n  G( K7 y! [dx(5)=x(6);
dx(6)=(400-(csp1*delta_sp11+ksp1*delta_sp1)*cosd(beita)*r1+T_in)/I1;
- n* X/ b$ j7 P! n# {dx(7)=x(8);
dx(8)=((-300*cosd(fai_sp1x)+200*cosd(fai_p1rx)-(cp1x*x(8))-kp1x*x(7)+(csp1*delta_sp11+ksp1*delta_sp1)*cosd(beita)*cosd(fai_sp1x)-(cp1r*delta_p11r+kp1r*delta_p1r))*cosd(beita)*cosd(fai_p1rx))/m2;                %第一个行星齿轮x方向
dx(9)=x(10);
dx(10)=((-300*cosd(fai_sp1y)+200*cosd(fai_p1ry)-(cp1y*x(10))-kp1y*x(9)+(csp1*delta_sp11+ksp1*delta_sp1)*cosd(beita)*cosd(fai_sp1y)+(cp1r*delta_p11r+kp1r*delta_p1r))*cosd(beita)*cosd(fai_p1ry))/m2;                 %第一个行星齿轮y方向
dx(11)=x(12);
dx(12)=(120-((csp1*delta_sp11+ksp1*delta_sp1)+(cp1r*delta_p11r+kp1r*delta_p1r))*cosd(beita)*r2)/I2;
dx(13)=x(14);
dx(14)=((-200*cosd(fai_p1rx)-(crx*x(14))-krx*x(13)+(cp1r*delta_p11r+kp1r*delta_p1r))*cosd(beita)*cosd(fai_p1rx))/m3;
! d/ M0 W4 ]8 y% E! {dx(15)=x(16);
dx(16)=((-200*cosd(fai_p1ry)-(crx*x(16))-krx*x(15)+(cp1r*delta_p11r+kp1r*delta_p1r))*cosd(beita)*cosd(fai_p1ry))/m3;
# Y6 H! f) T; _4 E' Bdx(17)=x(18);
dx(18)=(80-(cp1r*delta_p11r+kp1r*delta_p1r)*cosd(beita)*r3-T_out)/I3;

lkk

这么专业的问题还是求助淘宝吧