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

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

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

下面这个程序有点小问题，麻烦各位能否改一下，其中庞加莱图算出的数据总是不太多，而且画出的图有点问题
#include "stdio.h"
#include "math.h"
#include "conio.h"
#include<stdlib.h>
#define PI 3.1415926
#define Epsilon1.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,f,F,sita,clearance;               

static double sta,stb;
int n1;
static double k;
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=0.0;y=0.0;y=0.0;y=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,y);
      fprintf(y1y3,"%-16.8f%-16.8f\n",y,y);
      fprintf(y3y4,"%-16.8f%-16.8f\n",y,y);

      fprintf(y1t,"%-16.8f%-16.8f\n",tt,y);
      fprintf(y2t,"%-16.8f%-16.8f\n",tt,y);
      fprintf(y3t,"%-16.8f%-16.8f\n",tt,y);
      fprintf(y4t,"%-16.8f%-16.8f\n",tt,y);
    }
    //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,y);
   fprintf(pkly1y3,"%-16.8f%-16.8f\n",y,y);
   fprintf(pkly3y4,"%-16.8f%-16.8f\n",y,y);
   
       }
       NN++;
       hfa();
    }
   }
}
}

printf("%-16.8f\n",tt);
}

printf("Over!\n");
getch();
}

void cca()
{
int i;
sta=t;
for(i=1;i<=M;i++)sta=y;
sta=tt;
}
void hfa()
{
int i;
t=sta;
for(i=1;i<=M;i++)y=sta;
tt=sta;
}
void ccb()
{
int i;
stb=t;
for(i=1;i<=M;i++)stb=y;
stb=tt;
}
void hfb()
{
int i;
t=stb;
for(i=1;i<=M;i++)y=stb;
tt=stb;
}
void ccc()   /* only usefully to 2*PI/w poincare*/
{ int i;
stb=t;
for(i=1;i<=M;i++)stb=y;
stb=tt;
}
void hfc()
{ int i;t=stb;
for(i=1;i<=M;i++)y=stb;
tt=stb;
}
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;
    hh=hh/2.0;hh2=hh/2.0;
    hfb();
    temp=2*temp;
    for(l=1;l<=temp;l++)R_K();
    Delta=fabs(y-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=hh*f;
      k=y;
      y=k+1.0/2.0*k;
}
t=t+hh2;
equation();
for(i=1;i<=M;i++)
{
    k=hh*f;
      y=k+1.0/2.0*k;
}
equation();
for(i=1;i<=M;i++)
{
    k=hh*f;
      y=k+k;
}
t=t+hh2;
equation();
for(i=1;i<=M;i++)
    k=hh*f;
for(i=1;i<=M;i++)
    y=k+(k+2.0*k+2.0*k+k)/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,y);
   fprintf(pkly1y3,"%-16.8f%-16.8f\n",y,y);
   fprintf(pkly3y4,"%-16.8f%-16.8f\n",y,y);
   printf("%-16.8f\n",tt);
       }
       NN++;
       hfa();
    }
   }
}
}
}

void equation()
{

int i,j,k;

for(j=0;j<=1;j++)
      F=0.0;
for(k=1;k<=9;k++)
{
      sita=0.0;
      clearance=0.0;
}
for(i=1;i<=9;i++)
{
      sita=2*PI*(i-1)/9+w_c*t;
      clearance=y*sin(sita)+y*cos(sita)-game;
      if(clearance<=0)
          clearance=0;
         
      F=F+kb*pow(clearance,1.5)*cos(sita);
      F=F+kb*pow(clearance,1.5)*sin(sita);
}

f=y;
f=-(c/m)*y+W/m-F/m;
f=y;
f=-(c/m)*y-F/m;

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

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

其中图形如下

佩刀 发表于 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=2*PI*(i-1)/9+w_c*t;
      clearance=y*sin(sita)+y*cos(sita)-game;
      if(clearance<=0)
          clearance=0;
         
      F=F+kb*pow(clearance,1.5)*cos(sita);
      F=F+kb*pow(clearance,1.5)*sin(sita);


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

查看完整版本: 关于滚动轴承-转子系统非线性动力学程序方面的一些问题