刚才上传的附件，不方便看。下面是用C编写的runge kutta求解转子振动方程
不好意思，刚才上传附件，不方便大家看。下面是我用C编的求解轴承转子振动的runge-kutta法，大家来看一下。,为什么我求解了1000个周期，其转子中心位移越来越大呢，是不是结果发散了？请大家多批评，指点。谢谢！#include "stdlib.h"
#include "math.h"
void rkt2(t,h,y,n,eps,f)
void(*f)();
int n;
double t,h,eps,y[];
{ int m,i,j,k;
double hh,p,dt,x,tt,q,a,*g,*b,*c,*d,*e;
g=malloc(n*sizeof(double));
b=malloc(n*sizeof(double));
c=malloc(n*sizeof(double));
d=malloc(n*sizeof(double));
e=malloc(n*sizeof(double));
hh=h;
m=1;
p=1+eps;
x=t;
for(i=0; i<=n-1; i++)
c=y;
while(p>=eps)
{
a=hh/2.0;
a=a;
a=hh;a=hh;
for(i=0; i<=n-1;i++)
{
g=y;
y=c;
}
dt=h/m;
t=x;
for(j=0;j<=m-1; j++)
{
(*f)(t,y,n,d);
for(i=0;i<=n-1;i++)
{
b=y;
e=y;
}
for(k=0;k<=2;k++)
{
for(i=0;i<=n-1;i++)
{
y=e+a*d;
b=b+a*d/3.0;
}
tt=t+a;
(*f)(tt,y,n,d);
}
for(i=0; i<=n-1; i++)
y=b+hh*d/6.0;
t=t+dt;
}
p=0.0;
for(i=0; i<=n-1;i++)
{
q=fabs(y-g);
if(q>p) p=q;
}
hh=hh/2.0;
m=m+m;
}
free(g);
free(b);
free(c);
free(d);
free(e);
return;
}
**********************************************************************************************************
#include "stdio.h"
#include "math.h"
#include "stdlib.h"
#define PI 3.1415
#define g 9.8
int num_ball, n_rotor;//num_ball 滚珠的数量，n_rotor 转子转速
double damp,e,kn,w_rotor,w_c,F_out,m,gama,BN,T_vc,dt;
//damp 阻尼；e-偏心；kn-刚度；w_rotor 转子角速度；
//w_c 内圈转速；F_out-外力；m-转子质量；gama-游隙；
//BN,轴承参数，用以表示变刚度频率和转子转速的关系
//T_vc 变刚度的周期，60/(n_rotor*BN);%
//dt 周期内积分步长，一般一个周期取200个点
main()
{
int i,j;
void rkt2f(double,double [],int, double[]);
double t,h,eps,y;
damp=200.0;
e=0.0;
kn=7.055e9;
n_rotor=10000;
num_ball = 9;
F_out=6.0;
m=1.0;
gama=2e-5;
BN=3.6;
w_rotor=n_rotor*PI/30.0; //转速转换为角速度
w_c=0.4*w_rotor; //滚珠公转角速度
T_vc=50.0/(3.0*n_rotor); //变刚度周期T_vc=60/(n_rotor*3.6);
dt=1.0/(12*n_rotor); //一个周期取200个点0.005*T_vc=1/(12*n_rotor)
t=0.0;
h=dt;
printf("dt=",T_vc);
eps=1e-7;
y=0.0;y=1e-8;y=0.0;y=1e-8; //initial
printf("\n");
printf("t=%7.3f y=%e y=%e y=%e y=%e\n",t,y,y,y,y);
{
for(i=1;i<=200000;i++)//计算1000个周期
{ rkt2(t,h,y,4,eps,rkt2f);
// mrsn(t,h,y,4,eps,rkt2f);
t=t+h;
}
printf("t=%7.3f y=%e y=%e y=%e y=%e\n", t,y,y,y,y);
}
}
void rkt2f(t,y,n,d)
int n;
double t,y[],d[];
{
int i,j,k;
double *f;
double *sita; //滚珠的角度
double *clearance; //滚珠与内圈的间隙
t=t;
n=n;
f=(double*)malloc(2*sizeof(double));
sita=(double*)malloc(num_ball*sizeof(double));
clearance=(double*)malloc(num_ball*sizeof(double));
for(j=0;j<=1;j++)
f=0.0;
for(k=0;k<num_ball;k++)
{
sita=0.0;
clearance=0.0;
}
for(i=0;i<num_ball;i++)
{
sita=2*PI*(i-1)/num_ball+w_c*t;
clearance=y*sin(sita)+y*cos(sita)-gama;
if(clearance<=0)
clearance=0;
f=f+kn*pow(clearance,1.5)*cos(sita);
f=f+kn*pow(clearance,1.5)*sin(sita);
}
d=y;
d=e*w_rotor*w_rotor*cos(w_rotor*t)+g+F_out-damp*y-f-f;
d=y;
d=e*w_rotor*w_rotor*sin(w_rotor*t)-damp*y-f-f;
free(f);
free(sita);
free(clearance);
return;
} 最后一个for循环中，sita为sita；clearance应为clearance，不知道怎么回事，显示的不正确，望谅解！ 最有一个for 循环中,sita 应为sita(sita数组的第i 个元素) clearance应为clearance(clearance数组的第i个元素) 谢谢分享
页:
[1]