关于滚动轴承-转子系统非线性动力学程序方面的一些问题

佩刀 · 发表于 2017-12-14 16:16

马上注册，结交更多好友，享用更多功能，让你轻松玩转社区。

您需要登录才可以下载或查看，没有账号？我要加入

x

下面这个程序有点小问题，麻烦各位能否改一下，其中庞加莱图算出的数据总是不太多，而且画出的图有点问题

#include "stdio.h"
#include "math.h"
#include "conio.h"
#include<stdlib.h>
#define PI 3.1415926
#define Epsilon 1.0e-7
#define Epsilon1 1.0e-6
#define g 9.8
double TimeSup,xtTimeInf,pklTimeStart;
double M;
int NN;/* Di NN Ge Poincare JieMian*/
double h,w;
double G,e,b,w0,bc,wc,p;
double m,c,q,k1,kc,f1,u,s,kn;
double ss,W,kb,game,n_rotor,w_rotor,w_c,omega;
double hh,hh2;
double t,tt,phi;
static double y[10],f[10],F[2],sita[10],clearance[10];
static double sta[10],stb[10];
int n1;
static double k[10][5];
void cca(); void hfa(); void ccb();void hfb();
void ccc(); void hfc();
void R_K(); void chR_K(); void equation();
void pkl2pi();
main()
{
FILE *y1y2,*y1y3,*y3y4;
FILE *y1t,*y2t,*y3t,*y4t;
FILE *pkly1y2,*pkly1y3,*pkly3y4;
int i;
printf("Please wait... ...\n");
y1y2=fopen("3 y1y2.txt","w");
y1y3=fopen("3 y1y3.txt","w");
y3y4=fopen("3 y3y4.txt","w");
y1t=fopen("3 y1t.txt","w");
y2t=fopen("3 y2t.txt","w");
y3t=fopen("3 y3t.txt","w");
y4t=fopen("3 y4t.txt","w");
pkly1y2=fopen("3 pkly1y2.txt","w");
pkly1y3=fopen("3 pkly1y3.txt","w");
pkly3y4=fopen("3 pkly3y4.txt","w");
/***************** parameter start **************/
M=4;h=PI/512;
m=0.9;c=300.0;kb=7.055e9;game=20.0e-6;W=20.0;
n_rotor=21500;
w_rotor=n_rotor*PI/30.0; //×aËù×a»»Îa½ÇËù¶è
w0=sqrt(kb*sqrt(game)/(m));
w_c=0.4*w_rotor;
e=c/(m*w0);
ss=kb*pow(game,1.5);
omega=w_rotor/w0;
y[1]=0.0;y[2]=0.0;y[3]=0.0;y[4]=0.0;
TimeSup=1100.0;xtTimeInf=TimeSup-100;
pklTimeStart=0.0;
/***************** parameter end ****************/
hh=h;hh2=hh/2.0;
tt=0.0;t=0.0;NN=1;
equation();
cca();
while(tt<=TimeSup)
{
chR_K();
if(tt>xtTimeInf)
{
fprintf(y1y2,"%-16.8f%-16.8f\n",y[1],y[2]);
fprintf(y1y3,"%-16.8f%-16.8f\n",y[1],y[3]);
fprintf(y3y4,"%-16.8f%-16.8f\n",y[3],y[4]);
fprintf(y1t,"%-16.8f%-16.8f\n",tt,y[1]);
fprintf(y2t,"%-16.8f%-16.8f\n",tt,y[2]);
fprintf(y3t,"%-16.8f%-16.8f\n",tt,y[3]);
fprintf(y4t,"%-16.8f%-16.8f\n",tt,y[4]);
}
//if(tt>xtTimeInf)
if(tt-NN*2*PI<0.0&&-1*tt+NN*2*PI<=h)
{
cca();
{ int temp1=0;double temp2=1.0;
while(temp1!=1.0)
{ chR_K();
if(tt<NN*2*PI)ccc();
if(tt>NN*2*PI)
{
while(temp2>Epsilon)
{if(tt<NN*2*PI)
ccc();
if(tt>NN*2*PI)
{ hfc();
hh=hh/10.0;hh2=hh/2.0;
}
R_K();
temp2=fabs(tt-NN*2*PI);
}
hh=h;hh2=hh/2.0;temp1=1.0;
if(tt>xtTimeInf)
{
fprintf(pkly1y2,"%-16.8f%-16.8f\n",y[1],y[2]);
fprintf(pkly1y3,"%-16.8f%-16.8f\n",y[1],y[3]);
fprintf(pkly3y4,"%-16.8f%-16.8f\n",y[3],y[4]);
}
NN++;
hfa();
}
}
}
}
printf("%-16.8f\n",tt);
}
printf("Over!\n");
getch();
}
void cca()
{
int i;
sta[0]=t;
for(i=1;i<=M;i++)sta[i]=y[i];
sta[5]=tt;
}
void hfa()
{
int i;
t=sta[0];
for(i=1;i<=M;i++)y[i]=sta[i];
tt=sta[5];
}
void ccb()
{
int i;
stb[0]=t;
for(i=1;i<=M;i++)stb[i]=y[i];
stb[5]=tt;
}
void hfb()
{
int i;
t=stb[0];
for(i=1;i<=M;i++)y[i]=stb[i];
tt=stb[5];
}
void ccc() /* only usefully to 2*PI/w poincare*/
{ int i;
stb[0]=t;
for(i=1;i<=M;i++)stb[i]=y[i];
stb[5]=tt;
}
void hfc()
{ int i;t=stb[0];
for(i=1;i<=M;i++)y[i]=stb[i];
tt=stb[5];
}
void chR_K()
{
long l,temp=1;
double Delta=1.0;
double tempy2;
hh=h;hh2=hh/2.0;
ccb();
R_K();
while(Delta>Epsilon1)
{
tempy2=y[2];
hh=hh/2.0;hh2=hh/2.0;
hfb();
temp=2*temp;
for(l=1;l<=temp;l++)R_K();
Delta=fabs(y[2]-tempy2)/15.0;
}
hh=h;
hh2=hh/2.0;
}
void R_K()
{
int i;
tt=tt+hh;
equation();
for(i=1;i<=M;i++)
{
k[i][1]=hh*f[i];
k[i][0]=y[i];
y[i]=k[i][0]+1.0/2.0*k[i][1];
}
t=t+hh2;
equation();
for(i=1;i<=M;i++)
{
k[i][2]=hh*f[i];
y[i]=k[i][0]+1.0/2.0*k[i][2];
}
equation();
for(i=1;i<=M;i++)
{
k[i][3]=hh*f[i];
y[i]=k[i][0]+k[i][3];
}
t=t+hh2;
equation();
for(i=1;i<=M;i++)
k[i][4]=hh*f[i];
for(i=1;i<=M;i++)
y[i]=k[i][0]+(k[i][1]+2.0*k[i][2]+2.0*k[i][3]+k[i][4])/6.0;
}
void pkl2pi(FILE *pkly1y2,FILE *pkly1y3,FILE *pkly3y4)
{
if(tt-NN*2*PI<0.0&&-1*tt+NN*2*PI<=h)
{
cca();
{ int temp1=0;double temp2=1.0;
while(temp1!=1.0)
{ chR_K();
if(tt<NN*2*PI)ccc();
if(tt>NN*2*PI)
{
while(temp2>Epsilon)
{if(tt<NN*2*PI)
ccc();
if(tt>NN*2*PI)
{ hfc();
hh=hh/10.0;hh2=hh/2.0;
}
R_K();
temp2=fabs(tt-NN*2*PI);
}
hh=h;hh2=hh/2.0;temp1=1.0;
if(tt>pklTimeStart)
{
fprintf(pkly1y2,"%-16.8f%-16.8f\n",y[1],y[2]);
fprintf(pkly1y3,"%-16.8f%-16.8f\n",y[1],y[3]);
fprintf(pkly3y4,"%-16.8f%-16.8f\n",y[3],y[4]);
printf("%-16.8f\n",tt);
}
NN++;
hfa();
}
}
}
}
}
void equation()
{
int i,j,k;
for(j=0;j<=1;j++)
F[j]=0.0;
for(k=1;k<=9;k++)
{
sita[k]=0.0;
clearance[k]=0.0;
}
for(i=1;i<=9;i++)
{
sita[i]=2*PI*(i-1)/9+w_c*t;
clearance[i]=y[3]*sin(sita[i])+y[1]*cos(sita[i])-game;
if(clearance[i]<=0)
clearance[i]=0;
F[0]=F[0]+kb*pow(clearance[i],1.5)*cos(sita[i]);
F[1]=F[1]+kb*pow(clearance[i],1.5)*sin(sita[i]);
}
f[1]=y[2];
f[2]=-(c/m)*y[2]+W/m-F[0]/m;
f[3]=y[4];
f[4]=-(c/m)*y[4]-F[1]/m;

复制代码

佩刀 · 发表于 2017-12-14 16:18

本帖最后由佩刀于 2017-12-14 16:26 编辑

其中图形如下 QQ截图20171214162444.png

佩刀 · 发表于 2017-12-14 16:28

程序后面少了一个“｝”

CCTony · 发表于 2018-1-18 16:18

楼主，你好，我对你这个庞加莱图挺感兴趣的，能指导一下你的这个程序参数定义部分的含义吗？谢谢！

CCTony · 发表于 2018-1-18 16:21

M=4;h=PI/512;
m=0.9;c=300.0;kb=7.055e9;game=20.0e-6;W=20.0;

xusha889 · 发表于 2018-7-11 13:07

for(i=1;i<=9;i++)
  {
   sita[i]=2*PI*(i-1)/9+w_c*t;
   clearance[i]=y[3]*sin(sita[i])+y[1]*cos(sita[i])-game;
   if(clearance[i]<=0)
      clearance[i]=0;

   F[0]=F[0]+kb*pow(clearance[i],1.5)*cos(sita[i]);
   F[1]=F[1]+kb*pow(clearance[i],1.5)*sin(sita[i]);

我觉得这段程序与理论实际不符。

账号		自动登录	找回密码
密码			我要加入

[稳定性与分岔] 关于滚动轴承-转子系统非线性动力学程序方面的一些问题

马上注册，结交更多好友，享用更多功能，让你轻松玩转社区。