fei_udf.c

叉排管阵的流弹不稳定2维模型 FLUENT软件的UDF的C源代码 修正的隐式欧拉算法

/*C Code for 2 D Fluid elastic Instability Prediction of Triangle Tube Bundles*/
/*****************************************************************************
   2 degree Fluid elastic Instability Prediction of Triangle Tube Bundles
   modified implicit Euler method
predictor: y[p]=y[n]+f(time[n],y[n])*dtime
corrector: y[n+1]=y[n]+[f(time[n],y[n])+f(time[n+1],y[p])*dtime/2
*****************************************************************************/
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#define theta M_PI/3.0        /* angle at the seperation point */
#define den  1.225            /* air density  kg/m^3 */
#define E 201e9               /* Younge's Modulus of Elasticity of Tube, Pa */
#define IXX 3.22e 10         
/* Geometrical Moment of Inertia of Tube, m^4, i.e. for a round section, it is given by PI*d^4/64,                                in this case, d is 9.0mm */
#define Rho 4.204             /* Mass of Unit Length of Tube, kg/m */
#define L 0.268               /* L is the tube rod length, which is 268mm in this case */
#define Delta 0.003*2         /* Damping of Tube, which can be meassured by test */
#define K E*IXX/L/L/L         
/* Stiffness of Tube, kg/ms^2, i.e. for a cantilever, it is given by E*I/L^3.*/
#define Wr 2*M_PI*4.5 
/* First Order Natural Frequency of Tube, rad/s, i.e. for a cantilever with a concentrated mass at its free end, it can be computed by sqrt(K/(M+0.2357ml)), where M is the concentrated mass and ml is mass of cantilever. if f is given, then Wr is given by 2*PI*f  */
#define Am 0.020              /* the maximum possible amplitude of tube, m */
#define P  0.0892             /* tube pitch m */
#define H  0.077249           /* longitude space m */
#define D  0.0383             /* diameter of tube m */ 
#define Ug 1.5                /* air flow velocity across the minimum gap m/s */
/* the x coordinates of tube */
staticfloat CGC[28]={0.0,0.0892,0.1784,0.2676,0.3568,0.446,0.0446,0.1338,0.223,0.3122,0.4014,0.0,0.0892,0.1784,0.2676,0.3568,0.446,0.0446,0.1338,0.223,0.3122,0.4014,0.0,0.0892,0.1784,0.2676,0.3568,0.446,};
static float Vel_pre[28];
static float V_pre[28];
static float X_pre[28];
static float F_pre[28];
/* Modified Implicit Euler Algorithm */ 
int MIEA(X,V,F,dt,k) 
float X,V,F,dt; 
int k;
{ 
   int i; 
   float V_next,X_next,dF;   	 
/* predictor */
   	V_next=V+dt*(F V*Delta X*K)/Rho; /* f1=(fx x1*Delta/Rho x2*K/Rho) */
	 X_next=X+dt*V; /* f2=x1 */
/* corrector */
   dF=F F_pre[k];
   for(i=0;i<3;i++){
    V_next=V+0.5*dt*((F V*Delta X*K)+(F+2.5*dF V_next*Delta X_next*K))/Rho; 
	 X_next=X+0.5*dt*(V+V_next); 
	  if(V_next>Wr*Am){
       printf("the amplitude is too big... the proceed is terminated!");
        return  1;
      }
   }
   V_pre[k]=V_next; X_pre[k]=X_next;F_pre[k]=F;
}
main()
{
   int i,j,rm,n=5;
   float time,dt;
   float *r,rs,ru,rv,s,am;  /* generating random float_point values between [0,1] */
   float CGX[28],CGV[28],CGF[28]; /* displacement, velocity and force of tube array */
   float kl12,kl13,kr12,kr13; /* coefficients of vibrating displacement */
   float E2,E3,E4,E5,E6;    /* angle coefficients */
   float NE,WN,CW,C,CE;  /* the four vibrating displacements neibouring to the centre tube */
   float Dis,Vel,DV,Q;
   char str[8];
   FILE *fp[28];
/* initializing time, displacement, velocity and force */
   time=0.0;dt=0.01;
     for(i=0;i<28;i++){CGX[i]=0.0;CGV[i]=0.0;CGF[i]=0.0;}
/* initializing the displacement and position of each tube */  
rs=65536;ru=2053.0;rv=13849.0;am=0.0001; 
/* am is the assumed initial tiny amplitude of tube */
      for(i=0;i<28;i++){
        s=2.; r=&s;  /* to obtain different random seed r by changing s value */              
        *r=ru*(*r)+rv;rm=(int)(*r/rs);*r=*r rm*rs;
          CGX[i]=am*cos(2*M_PI*(*r/rs));
           CGV[i]=am*sin(2*M_PI*(*r/rs));
            CGC[i]=CGC[i]+CGX[i];
          }
/* constant coefficient calculations */ 
  E2=sin(theta); /* E2 */
  E3=((2*theta*theta 4)*sin(theta)+4*theta*cos(theta))/(theta*theta); /* E3 */
  E4=((theta*theta*theta*theta 12*theta*theta+24)*sin(theta)
     +(4*theta*theta*theta 24*theta)*cos(theta))/theta/theta/theta/theta; /* E4 */
  E5=sqrt((H D*sin(theta))*(H D*sin(theta))+D*D*cos(theta)*cos(theta)/4)*sin(theta); /* E5 */    
  E6=D*theta*sin(theta); /* E6 */
  /* gap mass flow rate */
  Q=den*Ug*(P D); 
  /*open files*/
   for(i=0;i<28;i++){
     sprintf(str,"%d",i+1); strcat(str,".txt");
       if((fp[i]=fopen(str,"at"))==NULL){printf("Cannot open the file!"); break;} 
    }
/* time step begin */    
   for(j=0;j<10000;j++){
/* Calculating the force on a tube */ 
   for(i=0;i<28;i++){
     /* first row */                
     if(i<6){
       WN=0.0;NE=0.0;
       if(i==0){CW=0.0;C=CGX[1];CE=CGX[2];} 
        else if(i==5){CW=CGX[4];C=CGX[5];CE=0.0;}
         else {CW=CGX[i 1];C=CGX[i];CE=CGX[i+1];} 
     }
     /* second row */ 
     else if(i>5&&i<11){
       WN=CGX[i 6];NE=CGX[i 5];
       if(i==6){CW=0.0;C=CGX[6];CE=CGX[7];} 
        else if(i==10){CW=CGX[9];C=CGX[10];CE=0.0;}
         else {CW=CGX[i 1];C=CGX[i];CE=CGX[i+1];}                    
     }
     /* third row */ 
     else if(i>10&&i<17){
       WN=CGX[i 6];NE=CGX[i 5];
       if(i==11){WN=0.0;CW=0.0;C=CGX[11];CE=CGX[12];} 
        else if(i==16){NE=0.0;CW=CGX[15];C=CGX[16];CE=0.0;}
         else {CW=CGX[i 1];C=CGX[i];CE=CGX[i+1];}                    
     }
     /* fourth row */ 
     else if(i>16&&i<22){
       WN=CGX[i 6];NE=CGX[i 5];
       if(i==17){CW=0.0;C=CGX[17];CE=CGX[18];} 
        else if(i==21){CW=CGX[20];C=CGX[21];CE=0.0;}
         else {CW=CGX[i 1];C=CGX[i];CE=CGX[i+1];}                    
     }
     /* fifth row */ 
     else {
       WN=CGX[i 6];NE=CGX[i 5];
       if(i==22){WN=0.0;CW=0.0;C=CGX[22];CE=CGX[23];} 
        else if(i==27){NE=0.0;CW=CGX[26];C=CGX[27];CE=0.0;}
         else {CW=CGX[i 1];C=CGX[i];CE=CGX[i+1];} 
     } 
     kl12=(P D*cos(theta)+NE WN)/(P D*cos(theta)+C CW); /* kl1 2*/
      kl13=(P D*cos(theta)+NE WN)/(P D+C CW); /* kl1 3*/
       kr12=(P D*cos(theta)+NE WN)/(P D*cos(theta)+CE C); /* kr1 2*/
        kr13=(P D*cos(theta)+NE WN)/(P D+CE C); /* kr1 3*/
     Dis=0.5*Ug*(P D)*D/(P D*cos(theta)+NE WN)/(P D*cos(theta)+NE WN);
       Dis=Dis*(E2*(kl13*kl13 kr13*kr13)+E3*((kl12 kl13)*kl13 (kr12 kr13)*kr13)
          +E4*((kl12 kl13)*(kl12 kl13) (kr12 kr13)*(kr12 kr13)));
     Vel=(E5+E6/3)*(kl12 kr12)+2*E6*(kl13 kr13)/3;
       Vel=0.5*Vel*D/(P D*cos(theta)+NE WN);
     DV=(Vel Vel_pre[i])/dt; Vel_pre[i]=Vel;
       CGF[i]= (Dis+DV)*Q; /* fluid force */
    /* computing the displacement and velocity of tube */
         rm=MIEA(CGX[i],CGV[i],CGF[i],dt,i);
          if(rm== 1)break;
    }
    for(i=0;i<28;i++){CGV[i]=V_pre[i];CGX[i]=X_pre[i];CGF[i]=F_pre[i];}
    /* write data file */
    if(j==n){
       for(i=0;i<28;i++)fprintf(fp[i], "%f  %f  %f  %f \n", time, V_pre[i],X_pre[i],F_pre[i]);
       n=n+2; }
     time=time+dt; }
    for(i=0;i<28;i++)fclose(fp[i]); }