关于用fluent UDF 编写风力机致动线模型



  • 各位老师好,我最近按照中科院王胜军博士的博士论文编写了一个风力机致动线模型的UDF,但一计算要么发散,要么浮点溢出。
    我的致动线的程序,我是把叶片分成40个叶素,先确每个叶素的中心点,风轮初始位置及运行的每一时刻,叶素中心的坐标好确定,但这个坐标不一定的网格中心的坐标,此时需要找到与之最近的网格作为线上的网格,并确定其中心坐标及相应网格中心的速度,然后返回去用BEM求解相应网格点的体积力,再用高斯分布分散到此网格周围,现在的问题的与叶素中心最近的那个网格怎么找呢?我用的一个遍历,可能不太对,还有我那个体积力源项可能也不太对,并行计算时,一开始就浮点溢出,串行时一开始就发散,我原来以为可能计算量大浮点溢出,但为什么串行也不行呢?所以我后来想着是我的UDDF出了问题,但是,我实在找不出问题出在哪里了。
    可能这个地方处理的不好
    efb3922e-d6ad-46e4-b4c1-0177ac9b5b09-image.png
    麻烦各位老师帮我看看,有没有更好一些的程序实现。
    程序流程图:
    e054e2b0-2390-4a2d-b283-14a087348e66-image.png
    这是部分程序代码,全部的好像传不上去,比较多。

    for(i=0;i<40;i++)             
        {
        rad[i]=(radius[i]+radius[i+1])/2;    /* 取分段点的中点为叶素*/
        drad[i]=rad[i+1]-rad[i];      /*叶素长度*/
        x01[i]=0;
        x02[i]=0;
        x03[i]=0;
        y01[i]=rad[i]*cos(thta0);
        y02[i]=rad[i]*cos(thta1);
        y03[i]=rad[i]*cos(thta2);
        z01[i]=rad[i]*sin(thta0);
        z02[i]=rad[i]*sin(thta1);
        z03[i]=rad[i]*sin(thta2);    
        }
      d1min=deta*2;
      d2min=deta*2; 
      d3min=deta*2; 
     
    for(i=0;i<40;i++)
       { 
        thread_loop_c(t,domain)      
        {          
         begin_c_loop_all(c,t)          
         {              
          C_CENTROID(m,c,t);
          x=m[0];
          y=m[1];
          z=m[2];
          d1=sqrt(pow(x-x01[i],2.0)+pow(y-y01[i],2.0)+pow(z-z01[i],2.0));
          d2=sqrt(pow(x-x02[i],2.0)+pow(y-y02[i],2.0)+pow(z-z02[i],2.0));
          d3=sqrt(pow(x-x03[i],2.0)+pow(y-y03[i],2.0)+pow(z-z03[i],2.0));
          if(d1<d1min) 
          {
           d1min=d1;
           x1[i]=x;
           y1[i]=y;
           z1[i]=z; 
          }
          if(d2<d1min) 
          {
           d2min=d2;
           x2[i]=x;
           y2[i]=y;
           z2[i]=z;
          }
          if(d3<d3min) 
          {
           d3min=d3;
           x3[i]=x;
           y3[i]=y;
           z3[i]=z;
          } 
          }
          end_c_loop(c,t)   
         }
      }
    


  • 部分程序01

    #include "udf.h" 
    #define pi 3.1415926 
    #define deta 0.171
    #define e 2.718281828459
    real dF1x[40],dF1y[40],dF1z[40],dF2x[40],dF2y[40],dF2z[40],dF3x[40],dF3y[40],dF3z[40],x1[40],x2[40],x3[40],y1[40],y2[40],y3[40],z1[40],z2[40],z3[40];
    DEFINE_INIT(my_init_function, domain) 
    { 
    Domain *d;
    cell_t c;    
    Thread *t; 
    real m[ND_ND];     /* loop over all cell threads in the domain   */  
    real flow_time=0;
    real visc,rho,tol,Numb,R,r_hub,U0,lambda,omg,thta0,thta1,thta2,x,y,z,radius[41],chord[40],twist[40],rad[40],drad[40];
    real sigma[40],a[40],b[40],newa,newb,phi[40],ft[40],fh[40],Ft[40],Fh[40],F[40],sinphi,cosphi,Re,aoe,Cl,Cd,Cn[40],Ct[40];
    real g1,g2,k,ac,aoa,dT[40],dM[40],u1[40],u2[40],u3[40],v1[40],v2[40],v3[40],w1[40],w2[40],w3[40];
    real  x01[40],y01[40],z01[40],x02[40],y02[40],z02[40],x03[40],y03[40],z03[40],d1min,d2min,d3min,d1,d2,d3;
    FILE *fpr,*fpc,*fpt;
    int i,j;
    
    visc=1.7894e-5; /* Kinematic viscosity of air (m^2/s) 运动粘度*/
    rho=1.225;     /* Density of air (kg/m^3) 密度*/
    tol=1.e-5;   /*Convergence tolerance for BE analysis 残差*/
    
    Numb=3;  /*叶片数*/
    R=7.4;   /*叶片半径*/
    r_hub=0.15;   /*轮毂半径*/
    
    U0=11;   /*额定风速*/
    lambda=6;  /*叶尖速比*/
    omg=lambda*U0/R;         /*弧度制转速*/ 
    
    thta0=omg*flow_time;
    thta1=omg*flow_time+pi*2/3;
    thta2=omg*flow_time+pi*4/3;
    
    fpr=fopen("F:\\mgl\\zhidongxian\\radius.txt","r");
    for(i=0;i<40;i++)	
       { 
       fscanf(fpr,"%f",&radius[i]);
        }								;
    fclose(fpr);
    fpc=fopen("F:\\mgl\\zhidongxian\\chord.txt","r");
    for(i=0;i<40;i++)	
       { 
       fscanf(fpc,"%f",&chord[i]);
        }								;
    fclose(fpc);
    fpt=fopen("F:\\mgl\\zhidongxian\\twist.txt","r");
    for(i=0;i<40;i++)	
       { 
       fscanf(fpt,"%f",&twist[i]);
        }								;
    fclose(fpt);
    radius[40]=7.4;
      
    for(i=0;i<40;i++)             
        {
        rad[i]=(radius[i]+radius[i+1])/2;    /* 取分段点的中点为叶素*/
        drad[i]=rad[i+1]-rad[i];      /*叶素长度*/
        x01[i]=0;
        x02[i]=0;
        x03[i]=0;
        y01[i]=rad[i]*cos(thta0);
        y02[i]=rad[i]*cos(thta1);
        y03[i]=rad[i]*cos(thta2);
        z01[i]=rad[i]*sin(thta0);
        z02[i]=rad[i]*sin(thta1);
        z03[i]=rad[i]*sin(thta2);    
        }
      d1min=deta*2;
      d2min=deta*2; 
      d3min=deta*2; 
     
    for(i=0;i<40;i++)
       { 
        thread_loop_c(t,domain)      
        {          
         begin_c_loop_all(c,t)          
         {              
          C_CENTROID(m,c,t);
          x=m[0];
          y=m[1];
          z=m[2];
          d1=sqrt(pow(x-x01[i],2.0)+pow(y-y01[i],2.0)+pow(z-z01[i],2.0));
          d2=sqrt(pow(x-x02[i],2.0)+pow(y-y02[i],2.0)+pow(z-z02[i],2.0));
          d3=sqrt(pow(x-x03[i],2.0)+pow(y-y03[i],2.0)+pow(z-z03[i],2.0));
          if(d1<d1min) 
          {
           d1min=d1;
           x1[i]=x;
           y1[i]=y;
           z1[i]=z; 
          }
          if(d2<d1min) 
          {
           d2min=d2;
           x2[i]=x;
           y2[i]=y;
           z2[i]=z;
          }
          if(d3<d3min) 
          {
           d3min=d3;
           x3[i]=x;
           y3[i]=y;
           z3[i]=z;
          } 
          }
          end_c_loop(c,t)   
         }
      }
        for(i=0;i<40;i++) 
       {        
           sigma[i]=Numb*chord[i]/(2*pi*radius[i]);  /*各截面实度*/    
           a[i]=0.3;b[i]=0.1;   /*给定a b 初始值*/
           newa=0.3;newb=0.1;
           for(j=0;j<4000;j++) 
           {
               phi[i]=atan(U0*(1-a[i])/(omg*radius[i]*(1+b[i])));
               ft[i]=Numb*(R-radius[i])/(2*radius[i]*sin(phi[i]));         /*叶尖损失系数参数*/
               fh[i]=Numb*(radius[i]- r_hub)/(2*r_hub*sin(phi[i]) ) ;   /*轮毂损失系数参数*/
               Ft[i]=2 *acos(exp(-ft[i]))/(pi);                     /*叶尖损失系数*/
               Fh[i]=2*acos(exp(-fh[i]))/(pi);                     /*轮毂损失系数*/
               F[i]=Ft[i]*Fh[i];           /*损失系数*/
               sinphi=sin(phi[i]);
               cosphi=cos(phi[i]);
    
               Re=chord[i]*rho* U0/visc;   /*雷诺数*/
               aoa=phi[i]*180/pi-twist[i];  /*角度制攻角*/
               if (i<2)                     /*Cylinder1*/
               {
               Cl=0; 
               Cd=0.5;                          /*获得Cl Cd*/
               }
               else if(i>=2&&i<4)
               {
               Cl=0;                              /*Cylinder2*/
               Cd=0.35;
                }
               else if(i>=4&&i<8)                    /*NACA443345E04*/
               {
                Cl=0.39176*pow(aoa,0)+0.16985*pow(aoa,1)-0.00948*pow(aoa,2)+2.05419E-4*pow(aoa,3)-1.98578E-6*pow(aoa,4)+7.04399E-9*pow(aoa,5);
                Cd=0.01299*pow(aoa,0)-0.00381*pow(aoa,1)+4.43552E-4*pow(aoa,2)-2.79172E-7*pow(aoa,3)-5.0702E-8*pow(aoa,4)+2.3948E-10*pow(aoa,5);
                }
               else if(i>=8&&i<12)                     /*NACA442545E04*/
               {
                Cl=0.33905*pow(aoa,0)+0.18374*pow(aoa,1)-0.01038*pow(aoa,2)+2.28113E-4*pow(aoa,3)-2.22776E-6*pow(aoa,4)+7.96503E-9*pow(aoa,5); 
                Cd=0.0163*pow(aoa,0)-0.00476*pow(aoa,1)+5.06172E-4*pow(aoa,2)-1.87549E-6*pow(aoa,3)-3.34699E-8*pow(aoa,4)+1.73246E-10*pow(aoa,5);
                }
               else if(i>=12&&i<16)                     /*NACA442345E04*/						
               {
                Cl=0.31081*pow(aoa,0)+0.18014*pow(aoa,1)-0.01003*pow(aoa,2)+2.18609E-4*pow(aoa,3)-2.12405E-6*pow(aoa,4)+7.56762E-9*pow(aoa,5); 
                Cd=0.01704*pow(aoa,0)-0.0048*pow(aoa,1)+5.04451E-4*pow(aoa,2)-1.79093E-6*pow(aoa,3)-3.4542E-8*pow(aoa,4)+1.77579E-10*pow(aoa,5); 
                }
               else if(i>=16&&i<23)                      /*NACA442265E04*/		
               {
                Cl=0.29245*pow(aoa,0)+0.17469*pow(aoa,1)-0.00961*pow(aoa,2)+2.07844E-4*pow(aoa,3)-2.01088E-6*pow(aoa,4)+7.14389E-9*pow(aoa,5); 
                Cd=0.01745*pow(aoa,0)-0.00471*pow(aoa,1)+5.00323E-4*pow(aoa,2)-1.71252E-6*pow(aoa,3)-3.52578E-8*pow(aoa,4)+1.80118E-10*pow(aoa,5);
                }
               else if(i>=23&&i<28)                        /*NACA442065E04*/	
               {
                Cl=0.29157*pow(aoa,0)+0.16896*pow(aoa,1)-0.00933*pow(aoa,2)+2.02874E-4*pow(aoa,3)-1.97245E-6*pow(aoa,4)+7.03532E-9*pow(aoa,5); 
                Cd=0.01862*pow(aoa,0)-0.00457*pow(aoa,1)+4.96343E-4*pow(aoa,2)-1.67345E-6*pow(aoa,3)-3.55333E-8*pow(aoa,4)+1.81392E-10*pow(aoa,5);
                }
               else if(i>=28&&i<35)                         /*NACA441865E04*/	
               {
                Cl=0.22608*pow(aoa,0)+0.16596*pow(aoa,1)-0.00896*pow(aoa,2)+1.92777E-4*pow(aoa,3)-1.86289E-6*pow(aoa,4)+6.61796E-9*pow(aoa,5); 
                Cd=0.01926*pow(aoa,0)-0.0047*pow(aoa,1)+4.9888E-4*pow(aoa,2)-1.66684E-6*pow(aoa,3)-3.59363E-8*pow(aoa,4)+1.83502E-10*pow(aoa,5);
                }
               else if(i>=35&&i<40)                             /*NACA441856E04*/
               {
                Cl=0.09669*pow(aoa,0)+0.14348*pow(aoa,1)-0.00739*pow(aoa,2)+1.57093E-4*pow(aoa,3)-1.51872E-6*pow(aoa,4)+5.41346E-9*pow(aoa,5); 
                Cd=0.02153*pow(aoa,0)-0.00387*pow(aoa,1)+4.72262E-4*pow(aoa,2)-1.32708E-6*pow(aoa,3)-3.82934E-8*pow(aoa,4)+1.91309E-10*pow(aoa,5);
               }
    
               Cn[i]=Cl*cosphi+Cd*sinphi ;  /*轴向力系数*/
               Ct[i]=Cl*sinphi-Cd*cosphi ;  /*切向力系数*/
    
               g1=sigma[i]*Cn[i]/(4*F[i]*sinphi*sinphi);
               g2=sigma[i]*Ct[i]/(4*F[i]*sinphi*cosphi);
              
               k=4*F[i]*sinphi*sinphi/(sigma[i]*Cn[i]);
               ac=0.38;
    
               if(a[i]<=0.38)       /*a 的修正*/
                 {    
                  newa=g1/(1+g1);                 
                  }
               else
                  { 
                  newa=0.5*(2+k*(1-2*ac)-sqrt(pow((k*(1-2*ac)+2),2)+4*(k*pow(ac,2) - 1))) ;
                  }
    
               newb=g2/(1-g2);
    
               if((fabs(newa-a[i])<=tol)&&(fabs(newb-b[i])<=tol))
                 { 
                  break;           
                  }
               else
                 {
               a[i]=newa;b[i]=newb;           
                  }
             }
           dT[i]=4*pi*rho*radius[i]*U0*U0*(1-a[i])*a[i]*F[i]*drad[i];   /*各叶素推力*/
           dM[i]=4*pi*rho*pow(radius[i],3)*U0*omg*(1-a[i])*b[i]*F[i]*drad[i]; /*各叶素转矩*/
           dF1x[i]=dT[i]; 
           dF2x[i]=dT[i];
           dF3x[i]=dT[i];
           dF1y[i]=0; 
           dF2y[i]=-(dM[i]/radius[i])*sin(thta1);
           dF3y[i]=-(dM[i]/radius[i])*sin(thta2);      
           dF1z[i]=dM[i]/radius[i];
           dF2z[i]=(dM[i]/radius[i])*cos(thta1);
           dF3z[i]=(dM[i]/radius[i])*cos(thta2);
        }                
    } 
    


  • 部分程序02

    DEFINE_SOURCE(ymom_source,c,t,dS,eqn)
    {
       real thta0,thta1,thta2,source,radius[41],chord[40],twist[40],rad[40],drad[40],visc,rho,tol,Numb,R,r_hub,U0,lambda,omg,flow_time;
       real epison,d,x,y,z;
       real eta3x1[40],eta3x2[40],eta3x3[40],eta3y1[40],eta3y2[40],eta3y3[40],eta3z1[40],eta3z2[40],eta3z3[40],dx1[40],dx2[40],dx3[40],dy1[40],dy2[40],dy3[40],dz1[40],dz2[40],dz3[40];
       real m[ND_ND];
       int i,n;
       FILE *fpr;
    
    visc=1.7894e-5; /* Kinematic viscosity of air (m^2/s) 运动粘滞系数*/
    rho=1.225;     /*Density of air (kg/m^3)密度*/ 
    tol=1.e-5;   /*Convergence tolerance for BE analysis 残差*/
    
    Numb=3;  /*叶片数*/
    R=7.4;   /*叶片半径*/
    r_hub=0.15;   /*轮毂半径*/
    
    U0=11;   /*额定风速*/
    lambda=6;  /*叶尖速比*/
    omg=lambda*U0/R;         /*弧度制转速*/ 
       
    flow_time=CURRENT_TIME;
    thta0=omg*flow_time;
    thta1=omg*flow_time+pi*2/3;
    thta2=omg*flow_time+pi*4/3;
    
    fpr=fopen("F:\\mgl\\zhidongxian\\radius.txt","r");
    for(i=0;i<40;i++)	
       { 
       fscanf(fpr,"%f",&radius[i]);
        }								;
    fclose(fpr);
    radius[40]=7.4;
    
    n=2;
    epison=n*deta;
    
          C_CENTROID(m,c,t);
          x=m[0];
          y=m[1];
          z=m[2]; 
            for(i=0;i<40;i++)
         { 
    	if(((x>=0)&&(x<deta))&&((pow(y,2.0)+pow(z,2.0))>=pow(radius[i],2))&&((pow(y,2.0)+pow(z,2.0))<=pow(radius[i+1],2)))
      	{
             y1[i]=y;
                 if((thta0==(pi/2+i*pi))) 
                      {
                       dy1[i]=y1[i];
                      } 
                 else 
                      {                                      
                      dy1[i]=((fabs(tan(-thta0)*y1[i]+z))/(sqrt(1+pow(tan(thta0),2.0))));
                      }
             eta3y1[i]=1/(pow(epison,3.)*pow(pi,3/2))*(exp(-pow(dy1[i]/epison,2.)));
    	 source=-eta3y1[i]*dF1y[i];   			
    	}       
            else if(((x>=0)&&(x<deta))&&((pow(y,2.0)+pow(z,2.0))>=pow(radius[i],2))&&((pow(y,2.0)+pow(z,2.0))<=pow(radius[i+1],2)))
      	{
             y2[i]=y;
                 if((thta1==(pi/2+i*pi))) 
                      {
                       dy2[i]=y2[i];
                      } 
                 else 
                      {                                      
                      dy2[i]=((fabs(tan(-thta1)*y2[i]+z))/(sqrt(1+pow(tan(thta1),2.0))));
                      }
             eta3y2[i]=1/(pow(epison,3.)*pow(pi,3/2))*(exp(-pow(dy2[i]/epison,2.)));
    	 source=-eta3y2[i]*dF2y[i];    
    	}
            else if(((x>=0)&&(x<deta))&&((pow(y,2.0)+pow(z,2.0))>=pow(radius[i],2))&&((pow(y,2.0)+pow(z,2.0))<=pow(radius[i+1],2)))
      	{
             y3[i]=y;
             if((thta2==(pi/2+i*pi))) 
                      {
                       dy3[i]=y3[i];
                      } 
                 else 
                      {                                      
                      dy3[i]=((fabs(tan(-thta2)*y3[i]+z))/(sqrt(1+pow(tan(thta2),2.0))));
                      }
             eta3y3[i]=1/(pow(epison,3.)*pow(pi,3/2))*(exp(-pow(dy3[i]/epison,2.)));
    	 source=-eta3y3[i]*dF3y[i];     
    	} 
          else
          {
           source=0;
           }
         } 
        return source;
    } 
    DEFINE_SOURCE(zmom_source,c,t,dS,eqn)
    {
       real thta0,thta1,thta2,source,radius[41],chord[40],twist[40],rad[40],drad[40],visc,rho,tol,Numb,R,r_hub,U0,lambda,omg,flow_time;
       real epison,d,x,y,z;
       real eta3x1[40],eta3x2[40],eta3x3[40],eta3y1[40],eta3y2[40],eta3y3[40],eta3z1[40],eta3z2[40],eta3z3[40],dx1[40],dx2[40],dx3[40],dy1[40],dy2[40],dy3[40],dz1[40],dz2[40],dz3[40];
       real m[ND_ND];
       int i,n;
       FILE *fpr;
    
    visc=1.7894e-5; /* Kinematic viscosity of air (m^2/s) 运动粘滞系数*/
    rho=1.225;     /*Density of air (kg/m^3)密度*/ 
    tol=1.e-5;   /*Convergence tolerance for BE analysis 残差*/
    
    Numb=3;  /*叶片数*/
    R=7.4;   /*叶片半径*/
    r_hub=0.15;   /*轮毂半径*/
    
    U0=11;   /*额定风速*/
    lambda=6;  /*叶尖速比*/
    omg=lambda*U0/R;         /*弧度制转速*/ 
       
    flow_time=CURRENT_TIME;
    thta0=omg*flow_time;
    thta1=omg*flow_time+pi*2/3;
    thta2=omg*flow_time+pi*4/3;
    
    fpr=fopen("F:\\mgl\\zhidongxian\\radius.txt","r");
    for(i=0;i<40;i++)	
       { 
       fscanf(fpr,"%f",&radius[i]);
        }								;
    fclose(fpr);
    radius[40]=7.4;
    
    n=2;
    epison=n*deta;
    
          C_CENTROID(m,c,t);
          x=m[0];
          y=m[1];
          z=m[2]; 
            for(i=0;i<40;i++)
         { 
    	if(((x>=0)&&(x<deta))&&((pow(y,2.0)+pow(z,2.0))>=pow(radius[i],2))&&((pow(y,2.0)+pow(z,2.0))<=pow(radius[i+1],2)))
      	{
    	 z1[i]=z;
                 if((thta0==(pi/2+i*pi))) 
                      {
                       dz1[i]=z1[i];
                      } 
                 else 
                      {                                      
                      dz1[i]=((fabs(tan(-thta0)*y+z1[i]))/(sqrt(1+pow(tan(thta0),2.0))));
                      }
             eta3z1[i]=1/(pow(epison,3.)*pow(pi,3/2))*(exp(-pow(dz1[i]/epison,2.)));
    	 source=-eta3z1[i]*dF1z[i];   				
    	}       
            else if(((x>=0)&&(x<deta))&&((pow(y,2.0)+pow(z,2.0))>=pow(radius[i],2))&&((pow(y,2.0)+pow(z,2.0))<=pow(radius[i+1],2)))
      	{
    	 z2[i]=z;
                 if((thta1==(pi/2+i*pi))) 
                      {
                       dz2[i]=z2[i];
                      } 
                 else 
                      {                                      
                      dz2[i]=((fabs(tan(-thta1)*y+z2[i]))/(sqrt(1+pow(tan(thta1),2.0))));
                      }
             eta3z2[i]=1/(pow(epison,3.)*pow(pi,3/2))*(exp(-pow(dz2[i]/epison,2.)));
    	 source=-eta3z2[i]*dF2z[i]; 
    	}
            else if(((x>=0)&&(x<deta))&&((pow(y,2.0)+pow(z,2.0))>=pow(radius[i],2))&&((pow(y,2.0)+pow(z,2.0))<=pow(radius[i+1],2)))
      	{
    	 z3[i]=z;
                 if((thta2==(pi/2+i*pi))) 
                      {
                       dz3[i]=z3[i];
                      } 
                 else 
                      {                                      
                      dz3[i]=((fabs(tan(-thta2)*y+z3[i]))/(sqrt(1+pow(tan(thta2),2.0))));
                      }
             eta3z3[i]=1/(pow(epison,3.0)*pow(pi,3/2))*(exp(-pow(dz3[i]/epison,2.0)));
    	 source=-eta3z3[i]*dF3z[i]; 
    	} 
          else
          {
           source=0;
           }
         } 
       return source;
     } 
    


  • 部分程序03

    DEFINE_ADJUST(my_adjust, domain) 
    {       
    Domain *d;
    cell_t c;    
    Thread *t; 
    real m[ND_ND];     /* loop over all cell threads in the domain   */  
    real flow_time=CURRENT_TIME;
    
    real visc,rho,tol,Numb,R,r_hub,U0,lambda,omg,thta0,thta1,thta2,x,y,z,radius[41],chord[40],twist[40],rad[40],drad[40];
    real sigma[40],a[40],b[40],newa,newb,phi[40],ft[40],fh[40],Ft[40],Fh[40],F[40],sinphi,cosphi,Re,aoe,Cl,Cd,Cn[40],Ct[40],g1,g2,k,ac;
    real dT[40],dM[40],u1[40],u2[40],u3[40],v1[40],v2[40],v3[40],w1[40],w2[40],w3[40];
    
    real a1[40],b1[40],newa1,newb1,a2[40],b2[40],newa2,newb2,a3[40],b3[40],newa3,newb3;
    real phi1[40],ft1[40],fh1[40],Ft1[40],Fh1[40],F1[40],phi2[40],ft2[40],fh2[40],Ft2[40],Fh2[40],F2[40];
    real phi3[40],ft3[40],fh3[40],Ft3[40],Fh3[40],F3[40],sinphi1,cosphi1,sinphi2,cosphi2,sinphi3,cosphi3;
    real aoa1,aoa2,aoa3,Cn1[40],Ct1[40],g11,g21,k1,dT1[40],dM1[40],Cn2[40],Ct2[40],g12,g22,k2,Re2,Re1,Re3;
    real dT2[40],dM2[40],Cn3[40],Ct3[40],g13,g23,k3,dT3[40],dM3[40],xiansudu1[40],xiansudu2[40],xiansudu3[40];
    real  x01[40],y01[40],z01[40],x02[40],y02[40],z02[40],x03[40],y03[40],z03[40],d1min,d2min,d3min,d1,d2,d3;
    
    FILE *fpr,*fpc,*fpt;
    int i,j;
    
    visc=1.7894e-5; /* Kinematic viscosity of air (m^2/s) 运动粘滞系数*/
    rho=1.225;     /*Density of air (kg/m^3)密度*/ 
    tol=1.e-5;   /*Convergence tolerance for BE analysis 残差*/
    
    Numb=3;  /*叶片数*/
    R=7.4;   /*叶片半径*/
    r_hub=0.15;   /*轮毂半径*/
    
    U0=11;   /*额定风速*/
    lambda=6;  /*叶尖速比*/
    omg=lambda*U0/R;         /*弧度制转速*/ 
    
    thta0=omg*flow_time;
    thta1=omg*flow_time+pi*2/3;
    thta2=omg*flow_time+pi*4/3;
    
    fpr=fopen("F:\\mgl\\zhidongxian\\radius.txt","r");
    for(i=0;i<40;i++)	
       { 
       fscanf(fpr,"%f",&radius[i]);
        }								;
    fclose(fpr);
    fpc=fopen("F:\\mgl\\zhidongxian\\chord.txt","r");
    for(i=0;i<40;i++)	
       { 
       fscanf(fpc,"%f",&chord[i]);
        }								
    fclose(fpc);
    fpt=fopen("F:\\mgl\\zhidongxian\\twist.txt","r");
    for(i=0;i<40;i++)	
       { 
       fscanf(fpt,"%f",&twist[i]);
        }								
    fclose(fpt);
    radius[40]=7.4;
    
    for(i=0;i<40;i++)             
        {
        rad[i]=(radius[i]+radius[i+1])/2;    /* 取分段点的中点为叶素*/
        drad[i]=rad[i+1]-rad[i];      /*叶素长度*/
        x01[i]=0;
        x02[i]=0;
        x03[i]=0;
        y01[i]=rad[i]*cos(thta0);
        y02[i]=rad[i]*cos(thta1);
        y03[i]=rad[i]*cos(thta2);
        z01[i]=rad[i]*sin(thta0);
        z02[i]=rad[i]*sin(thta1);
        z03[i]=rad[i]*sin(thta2);    
        }
      d1min=deta*2;
      d2min=deta*2; 
      d3min=deta*2; 
    
    for(i=0;i<40;i++)
       { 
        thread_loop_c(t,domain)      
        {          
         begin_c_loop_all(c,t)          
         {              
          C_CENTROID(m,c,t);
          x=m[0];
          y=m[1];
          z=m[2];
          d1=sqrt(pow(x-x01[i],2.0)+pow(y-y01[i],2.0)+pow(z-z01[i],2.0));
          d2=sqrt(pow(x-x02[i],2.0)+pow(y-y02[i],2.0)+pow(z-z02[i],2.0));
          d3=sqrt(pow(x-x03[i],2.0)+pow(y-y03[i],2.0)+pow(z-z03[i],2.0));
          if(d1<d1min) 
          {
           d1min=d1;
           x1[i]=x;
           y1[i]=y;
           z1[i]=z;
           u1[i]=C_U(c,t);
           v1[i]=C_V(c,t);
           w1[i]=C_W(c,t); 
          }
          if(d2<d2min) 
          {
           d2min=d2;
           x2[i]=x;
           y2[i]=y;
           z2[i]=z;
           u2[i]=C_U(c,t);
           v2[i]=C_V(c,t);
           w2[i]=C_W(c,t);
          }
          if(d3<d3min) 
          {
           d3min=d3;
           x3[i]=x;
           y3[i]=y;
           z3[i]=z;
           u2[i]=C_U(c,t);
           v2[i]=C_V(c,t);
           w2[i]=C_W(c,t);
          } 
          }
          end_c_loop(c,t)   
         } 
    }      				
        for(i=0;i<40;i++) 
      {
           sigma[i]=Numb*chord[i]/(2*pi*radius[i]);  /*各截面实度*/    
           a1[i]=0.3;b1[i]=0.1;   /*给定a b 初始值*/
           newa1=0.3;newb1=0.1;
           a2[i]=0.3;b2[i]=0.1;   /*给定a b 初始值*/
           newa2=0.3;newb2=0.1;
           a3[i]=0.3;b3[i]=0.1;   /*给定a b 初始值*/
           newa3=0.3;newb3=0.1;
           xiansudu1[i]=sqrt(v1[i]+w1[i]);
           xiansudu2[i]=sqrt(v2[i]+w2[i]);
           xiansudu3[i]=sqrt(v3[i]+w3[i]);
          for(j=0;j<4000;j++) 
           {
               phi1[i]=atan(u1[i]*(1-a1[i])/(xiansudu1[i]*(1+b1[i])));
               ft1[i]=Numb*(R-radius[i])/(2*radius[i]*sin(phi1[i]));/*叶尖损失系数参数*/                                                                    
               fh1[i]=Numb*(radius[i]- r_hub)/( 2*r_hub*sin(phi1[i]));   /*轮毂损失系数参数*/
               Ft1[i]=2 *acos(exp(-ft1[i]))/(pi);                     /*叶尖损失系数*/
               Fh1[i]=2*acos(exp(-fh1[i]))/(pi);                     /*轮毂损失系数*/  
               F1[i]=Ft1[i]*Fh1[i];           /*损失系数*/
    
               sinphi1=sin(phi1[i]);
               cosphi1=cos(phi1[i]);
    
               Re1=chord[i]*rho*u1[i]/visc;   /*雷诺数*/
               aoa1=phi1[i]*180/pi-twist[i];  /*角度制攻角*/
               if (i<2)                     /*Cylinder1*/
               {
               Cl=0; 
               Cd=0.5;                          /*获得Cl Cd*/
               }
               else if(i>=2&&i<4)
               {
               Cl=0;                              /*Cylinder2*/
               Cd=0.35;
                }
               else if(i>=4&&i<8)                    /*NACA443345E04*/
               {
                Cl=0.39176*pow(aoa1,0)+0.16985*pow(aoa1,1)-0.00948*pow(aoa1,2)+2.05419E-4*pow(aoa1,3)-1.98578E-6*pow(aoa1,4)+7.04399E-9*pow(aoa1,5);
                Cd=0.01299*pow(aoa1,0)-0.00381*pow(aoa1,1)+4.43552E-4*pow(aoa1,2)-2.79172E-7*pow(aoa1,3)-5.0702E-8*pow(aoa1,4)+2.3948E-10*pow(aoa1,5);
                }
               else if(i>=8&&i<12)                     /*NACA442545E04*/
               {
                Cl=0.33905*pow(aoa1,0)+0.18374*pow(aoa1,1)-0.01038*pow(aoa1,2)+2.28113E-4*pow(aoa1,3)-2.22776E-6*pow(aoa1,4)+7.96503E-9*pow(aoa1,5); 
                Cd=0.0163*pow(aoa1,0)-0.00476*pow(aoa1,1)+5.06172E-4*pow(aoa1,2)-1.87549E-6*pow(aoa1,3)-3.34699E-8*pow(aoa1,4)+1.73246E-10*pow(aoa1,5);
                }
               else if(i>=12&&i<16)                     /*NACA442345E04*/						
               {
                Cl=0.31081*pow(aoa1,0)+0.18014*pow(aoa1,1)-0.01003*pow(aoa1,2)+2.18609E-4*pow(aoa1,3)-2.12405E-6*pow(aoa1,4)+7.56762E-9*pow(aoa1,5); 
                Cd=0.01704*pow(aoa1,0)-0.0048*pow(aoa1,1)+5.04451E-4*pow(aoa1,2)-1.79093E-6*pow(aoa1,3)-3.4542E-8*pow(aoa1,4)+1.77579E-10*pow(aoa1,5); 
                }
               else if(i>=16&&i<23)                      /*NACA442265E04*/		
               {
                Cl=0.29245*pow(aoa1,0)+0.17469*pow(aoa1,1)-0.00961*pow(aoa1,2)+2.07844E-4*pow(aoa1,3)-2.01088E-6*pow(aoa1,4)+7.14389E-9*pow(aoa1,5); 
                Cd=0.01745*pow(aoa1,0)-0.00471*pow(aoa1,1)+5.00323E-4*pow(aoa1,2)-1.71252E-6*pow(aoa1,3)-3.52578E-8*pow(aoa1,4)+1.80118E-10*pow(aoa1,5);
                }
               else if(i>=23&&i<28)                        /*NACA442065E04*/	
               {
                Cl=0.29157*pow(aoa1,0)+0.16896*pow(aoa1,1)-0.00933*pow(aoa1,2)+2.02874E-4*pow(aoa1,3)-1.97245E-6*pow(aoa1,4)+7.03532E-9*pow(aoa1,5); 
                Cd=0.01862*pow(aoa1,0)-0.00457*pow(aoa1,1)+4.96343E-4*pow(aoa1,2)-1.67345E-6*pow(aoa1,3)-3.55333E-8*pow(aoa1,4)+1.81392E-10*pow(aoa1,5);
                }
               else if(i>=28&&i<35)                         /*NACA441865E04*/	
               {
                Cl=0.22608*pow(aoa1,0)+0.16596*pow(aoa1,1)-0.00896*pow(aoa1,2)+1.92777E-4*pow(aoa1,3)-1.86289E-6*pow(aoa1,4)+6.61796E-9*pow(aoa1,5); 
                Cd=0.01926*pow(aoa1,0)-0.0047*pow(aoa1,1)+4.9888E-4*pow(aoa1,2)-1.66684E-6*pow(aoa1,3)-3.59363E-8*pow(aoa1,4)+1.83502E-10*pow(aoa1,5);
                }
               else if(i>=35&&i<40)                             /*NACA441856E04*/
               {
                Cl=0.09669*pow(aoa1,0)+0.14348*pow(aoa1,1)-0.00739*pow(aoa1,2)+1.57093E-4*pow(aoa1,3)-1.51872E-6*pow(aoa1,4)+5.41346E-9*pow(aoa1,5); 
                Cd=0.02153*pow(aoa1,0)-0.00387*pow(aoa1,1)+4.72262E-4*pow(aoa1,2)-1.32708E-6*pow(aoa1,3)-3.82934E-8*pow(aoa1,4)+1.91309E-10*pow(aoa1,5);
               }
    
               Cn1[i]=Cl*cosphi1+Cd*sinphi1 ;   /*轴向力系数*/
               Ct1[i]=Cl*sinphi1-Cd*cosphi1 ;   /*切向力系数*/
    
               g11=sigma[i]*Cn1[i]/(4*F1[i]*sinphi1*sinphi1);
               g21=sigma[i]*Ct1[i]/(4*F1[i]*sinphi1*cosphi1);
              
               k1=4*F1[i]*sinphi1*sinphi1/(sigma[i]*Cn1[i]);
               ac=0.38;
    
               if(a1[i]<=0.38)       /*a 的修正*/
                 {    
                  newa1=g11/(1+g11);                 
                  }
               else
                  { 
                  newa1=0.5*(2+k1*(1-2*ac)-sqrt(pow((k1*(1-2*ac)+2),2)+4*(k1*pow(ac,2) - 1))) ;
                  }
    
               newb1=g21/(1-g21);
    
               if((fabs(newa1-a1[i])<=tol)&&(fabs(newb1-b1[i])<=tol))
                 { 
                  break;           
                  }
               else
                 {
                 a1[i]=newa1;b1[i]=newb1;           
                  }
            }
           dT1[i]=4*pi*rho*radius[i]*u1[i]*u1[i]*(1-a1[i])*a1[i]*F1[i]*drad[i];   /*各叶素推力*/
           dM1[i]=4*pi*rho*pow(radius[i],3)*u1[i]*xiansudu1[i]/radius[i]*(1-a1[i])*b1[i]*F1[i]*drad[i]; /*各叶素转矩*/
         
          for(j=0;j<4000;j++) 
           {
               phi2[i]=atan(u2[i]*(1-a2[i])/(xiansudu2[i]*(1+b2[i])));
               ft2[i]=Numb*(R-radius[i])/(2*radius[i]*sin(phi2[i]));/*叶尖损失系数参数*/                                                                    
               fh2[i]=Numb*(radius[i]-r_hub)/( 2*r_hub*sin(phi2[i])) ;   /*轮毂损失系数参数*/
               Ft2[i]=2 *acos(exp(-ft2[i]))/(pi);                     /*叶尖损失系数*/
               Fh2[i]=2*acos(exp(-fh2[i]))/(pi);                     /*轮毂损失系数*/  
               F2[i]=Ft2[i]*Fh2[i];           /*损失系数*/
    
               sinphi2=sin(phi2[i]);
               cosphi2=cos(phi2[i]);
    
               Re2=chord[i]*rho*u2[i]/visc;   /*雷诺数*/
               aoa2=phi2[i]*180/pi-twist[i];  /*角度制攻角*/
               if (i<2)                     /*Cylinder1*/
               {
               Cl=0; 
               Cd=0.5;                          /*获得Cl Cd*/
               }
               else if(i>=2&&i<4)
               {
               Cl=0;                              /*Cylinder2*/
               Cd=0.35;
                }
               else if(i>=4&&i<8)                    /*NACA443345E04*/
               {
                Cl=0.39176*pow(aoa2,0)+0.16985*pow(aoa2,1)-0.00948*pow(aoa2,2)+2.05419E-4*pow(aoa2,3)-1.98578E-6*pow(aoa2,4)+7.04399E-9*pow(aoa2,5);
                Cd=0.01299*pow(aoa2,0)-0.00381*pow(aoa2,1)+4.43552E-4*pow(aoa2,2)-2.79172E-7*pow(aoa2,3)-5.0702E-8*pow(aoa2,4)+2.3948E-10*pow(aoa2,5);
                }
               else if(i>=8&&i<12)                     /*NACA442545E04*/
               {
                Cl=0.33905*pow(aoa2,0)+0.18374*pow(aoa2,1)-0.01038*pow(aoa2,2)+2.28113E-4*pow(aoa2,3)-2.22776E-6*pow(aoa2,4)+7.96503E-9*pow(aoa2,5); 
                Cd=0.0163*pow(aoa2,0)-0.00476*pow(aoa2,1)+5.06172E-4*pow(aoa2,2)-1.87549E-6*pow(aoa2,3)-3.34699E-8*pow(aoa2,4)+1.73246E-10*pow(aoa2,5);
                }
               else if(i>=12&&i<16)                     /*NACA442345E04*/						
               {
                Cl=0.31081*pow(aoa2,0)+0.18014*pow(aoa2,1)-0.01003*pow(aoa2,2)+2.18609E-4*pow(aoa2,3)-2.12405E-6*pow(aoa2,4)+7.56762E-9*pow(aoa2,5); 
                Cd=0.01704*pow(aoa2,0)-0.0048*pow(aoa2,1)+5.04451E-4*pow(aoa2,2)-1.79093E-6*pow(aoa2,3)-3.4542E-8*pow(aoa2,4)+1.77579E-10*pow(aoa2,5); 
                }
               else if(i>=16&&i<23)                      /*NACA442265E04*/		
               {
                Cl=0.29245*pow(aoa2,0)+0.17469*pow(aoa2,1)-0.00961*pow(aoa2,2)+2.07844E-4*pow(aoa2,3)-2.01088E-6*pow(aoa2,4)+7.14389E-9*pow(aoa2,5); 
                Cd=0.01745*pow(aoa2,0)-0.00471*pow(aoa2,1)+5.00323E-4*pow(aoa2,2)-1.71252E-6*pow(aoa2,3)-3.52578E-8*pow(aoa2,4)+1.80118E-10*pow(aoa2,5);
                }
               else if(i>=23&&i<28)                        /*NACA442065E04*/	
               {
                Cl=0.29157*pow(aoa2,0)+0.16896*pow(aoa2,1)-0.00933*pow(aoa2,2)+2.02874E-4*pow(aoa2,3)-1.97245E-6*pow(aoa2,4)+7.03532E-9*pow(aoa2,5); 
                Cd=0.01862*pow(aoa2,0)-0.00457*pow(aoa2,1)+4.96343E-4*pow(aoa2,2)-1.67345E-6*pow(aoa2,3)-3.55333E-8*pow(aoa2,4)+1.81392E-10*pow(aoa2,5);
                }
               else if(i>=28&&i<35)                         /*NACA441865E04*/	
               {
                Cl=0.22608*pow(aoa2,0)+0.16596*pow(aoa2,1)-0.00896*pow(aoa2,2)+1.92777E-4*pow(aoa2,3)-1.86289E-6*pow(aoa2,4)+6.61796E-9*pow(aoa2,5); 
                Cd=0.01926*pow(aoa2,0)-0.0047*pow(aoa2,1)+4.9888E-4*pow(aoa2,2)-1.66684E-6*pow(aoa2,3)-3.59363E-8*pow(aoa2,4)+1.83502E-10*pow(aoa2,5);
                }
               else if(i>=35&&i<40)                             /*NACA441856E04*/
               {
                Cl=0.09669*pow(aoa2,0)+0.14348*pow(aoa2,1)-0.00739*pow(aoa2,2)+1.57093E-4*pow(aoa2,3)-1.51872E-6*pow(aoa2,4)+5.41346E-9*pow(aoa2,5); 
                Cd=0.02153*pow(aoa2,0)-0.00387*pow(aoa2,1)+4.72262E-4*pow(aoa2,2)-1.32708E-6*pow(aoa2,3)-3.82934E-8*pow(aoa2,4)+1.91309E-10*pow(aoa2,5);
               }
    
               Cn2[i]=Cl*cosphi2+Cd*sinphi2 ;   /*轴向力系数*/
               Ct2[i]=Cl*sinphi2-Cd*cosphi2 ;   /*切向力系数*/
    
               g12=sigma[i]*Cn1[i]/(4*F2[i]*sinphi2*sinphi2);
               g22=sigma[i]*Ct1[i]/(4*F2[i]*sinphi2*cosphi2);
              
               k2=4*F2[i]*sinphi2*sinphi2/(sigma[i]*Cn2[i]);
               ac=0.38;
    
               if(a2[i]<=0.38)       /*a 的修正*/
                 {    
                  newa1=g12/(1+g12);                 
                  }
               else
                  { 
                  newa1=0.5*(2+k2*(1-2*ac)-sqrt(pow((k2*(1-2*ac)+2),2)+4*(k2*pow(ac,2) - 1))) ;
                  }
    
               newb1=g22/(1-g22);
    
               if((fabs(newa2-a2[i])<=tol)&&(fabs(newb2-b2[i])<=tol))
                 { 
                  break;           
                  }
               else
                 {
               a2[i]=newa2;b2[i]=newb2;           
                 }
          }
           dT2[i]=4*pi*rho*radius[i]*u2[i]*u2[i]*(1-a2[i])*a2[i]*F2[i]*drad[i];   /*各叶素推力*/
           dM2[i]=4*pi*rho*pow(radius[i],3)*u2[i]*xiansudu2[i]/radius[i]*(1-a2[i])*b2[i]*F2[i]*drad[i]; /*各叶素转矩*/
    
        for(j=0;j<4000;j++) 
           {
               phi3[i]=atan(u3[i]*(1-a3[i])/(xiansudu3[i]*(1+b3[i])));
               ft3[i]=Numb*(R-radius[i])/(2*radius[i]*sin(phi3[i]));/*叶尖损失系数参数*/                                                                    
               fh3[i]=Numb*(radius[i]- r_hub)/( 2*r_hub*sin(phi3[i])) ;   /*轮毂损失系数参数*/
               Ft3[i]=2 *acos(exp(-ft3[i]))/(pi);                     /*叶尖损失系数*/
               Fh3[i]=2*acos(exp(-fh3[i]))/(pi);                     /*轮毂损失系数*/  
               F3[i]=Ft3[i]*Fh3[i];           /*损失系数*/
    
               sinphi3=sin(phi3[i]);
               cosphi3=cos(phi3[i]);
    
               Re3=chord[i]*rho*u3[i]/visc;   /*雷诺数*/
               aoa3=phi3[i]*180/pi-twist[i];  /*角度制攻角*/
               if (i<2)                     /*Cylinder1*/
               {
               Cl=0; 
               Cd=0.5;                          /*获得Cl Cd*/
               }
               else if(i>=2&&i<4)
               {
               Cl=0;                              /*Cylinder2*/
               Cd=0.35;
                }
               else if(i>=4&&i<8)                    /*NACA443345E04*/
               {
                Cl=0.39176*pow(aoa3,0)+0.16985*pow(aoa3,1)-0.00948*pow(aoa3,2)+2.05419E-4*pow(aoa3,3)-1.98578E-6*pow(aoa3,4)+7.04399E-9*pow(aoa3,5);
                Cd=0.01299*pow(aoa3,0)-0.00381*pow(aoa3,1)+4.43552E-4*pow(aoa3,2)-2.79172E-7*pow(aoa3,3)-5.0702E-8*pow(aoa3,4)+2.3948E-10*pow(aoa3,5);
                }
               else if(i>=8&&i<12)                     /*NACA442545E04*/
               {
                Cl=0.33905*pow(aoa3,0)+0.18374*pow(aoa3,1)-0.01038*pow(aoa3,2)+2.28113E-4*pow(aoa3,3)-2.22776E-6*pow(aoa3,4)+7.96503E-9*pow(aoa3,5); 
                Cd=0.0163*pow(aoa3,0)-0.00476*pow(aoa3,1)+5.06172E-4*pow(aoa3,2)-1.87549E-6*pow(aoa3,3)-3.34699E-8*pow(aoa3,4)+1.73246E-10*pow(aoa3,5);
                }
               else if(i>=12&&i<16)                     /*NACA442345E04*/						
               {
                Cl=0.31081*pow(aoa3,0)+0.18014*pow(aoa3,1)-0.01003*pow(aoa3,2)+2.18609E-4*pow(aoa3,3)-2.12405E-6*pow(aoa3,4)+7.56762E-9*pow(aoa3,5); 
                Cd=0.01704*pow(aoa3,0)-0.0048*pow(aoa3,1)+5.04451E-4*pow(aoa3,2)-1.79093E-6*pow(aoa3,3)-3.4542E-8*pow(aoa3,4)+1.77579E-10*pow(aoa3,5); 
                }
               else if(i>=16&&i<23)                      /*NACA442265E04*/		
               {
                Cl=0.29245*pow(aoa3,0)+0.17469*pow(aoa3,1)-0.00961*pow(aoa3,2)+2.07844E-4*pow(aoa3,3)-2.01088E-6*pow(aoa3,4)+7.14389E-9*pow(aoa3,5); 
                Cd=0.01745*pow(aoa3,0)-0.00471*pow(aoa3,1)+5.00323E-4*pow(aoa3,2)-1.71252E-6*pow(aoa3,3)-3.52578E-8*pow(aoa3,4)+1.80118E-10*pow(aoa3,5);
                }
               else if(i>=23&&i<28)                        /*NACA442065E04*/	
               {
                Cl=0.29157*pow(aoa3,0)+0.16896*pow(aoa3,1)-0.00933*pow(aoa3,2)+2.02874E-4*pow(aoa3,3)-1.97245E-6*pow(aoa3,4)+7.03532E-9*pow(aoa3,5); 
                Cd=0.01862*pow(aoa3,0)-0.00457*pow(aoa3,1)+4.96343E-4*pow(aoa3,2)-1.67345E-6*pow(aoa3,3)-3.55333E-8*pow(aoa3,4)+1.81392E-10*pow(aoa3,5);
                }
               else if(i>=28&&i<35)                         /*NACA441865E04*/	
               {
                Cl=0.22608*pow(aoa3,0)+0.16596*pow(aoa3,1)-0.00896*pow(aoa3,2)+1.92777E-4*pow(aoa3,3)-1.86289E-6*pow(aoa3,4)+6.61796E-9*pow(aoa3,5); 
                Cd=0.01926*pow(aoa3,0)-0.0047*pow(aoa3,1)+4.9888E-4*pow(aoa3,2)-1.66684E-6*pow(aoa3,3)-3.59363E-8*pow(aoa3,4)+1.83502E-10*pow(aoa3,5);
                }
               else if(i>=35&&i<40)                             /*NACA441856E04*/
               {
                Cl=0.09669*pow(aoa3,0)+0.14348*pow(aoa3,1)-0.00739*pow(aoa3,2)+1.57093E-4*pow(aoa3,3)-1.51872E-6*pow(aoa3,4)+5.41346E-9*pow(aoa3,5); 
                Cd=0.02153*pow(aoa3,0)-0.00387*pow(aoa3,1)+4.72262E-4*pow(aoa3,2)-1.32708E-6*pow(aoa3,3)-3.82934E-8*pow(aoa3,4)+1.91309E-10*pow(aoa3,5);
               }
    
               Cn3[i]=Cl*cosphi3+Cd*sinphi3 ;   /*轴向力系数*/
               Ct3[i]=Cl*sinphi3-Cd*cosphi3 ;   /*切向力系数*/
    
               g13=sigma[i]*Cn3[i]/(4*F2[i]*sinphi3*sinphi3);
               g23=sigma[i]*Ct3[i]/(4*F2[i]*sinphi3*cosphi3);
              
               k3=4*F2[i]*sinphi3*sinphi3/(sigma[i]*Cn3[i]);
               ac=0.38;
    
               if(a3[i]<=0.38)       /*a 的修正*/
                 {    
                  newa3=g13/(1+g13);                 
                  }
               else
                  { 
                  newa3=0.5*(2+k3*(1-2*ac)-sqrt(pow((k3*(1-2*ac)+2),2)+4*(k3*pow(ac,2) - 1))) ;
                  }
    
               newb3=g23/(1-g23);
    
               if((fabs(newa3-a3[i])<=tol)&&(fabs(newb3-b3[i])<=tol))
                 { 
                  break;           
                  }
               else
                 {
               a3[i]=newa3;b3[i]=newb3;           
                 }
            }
           dT3[i]=4*pi*rho*radius[i]*u3[i]*u3[i]*(1-a3[i])*a3[i]*F3[i]*drad[i];   /*各叶素推力*/
           dM3[i]=4*pi*rho*pow(radius[i],3)*u3[i]*xiansudu3[i]/radius[i]*(1-a3[i])*b3[i]*F3[i]*drad[i]; /*各叶素转矩*/
        }      
           dF1x[i]=dT1[i]; 
           dF2x[i]=dT2[i];
           dF3x[i]=dT3[i];
           dF1y[i]=-(dM1[i]/radius[i])*sin(thta0); 
           dF2y[i]=-(dM2[i]/radius[i])*sin(thta1);
           dF3y[i]=-(dM3[i]/radius[i])*sin(thta2);      
           dF1z[i]=(dM1[i]/radius[i])*cos(thta0);
           dF2z[i]=(dM2[i]/radius[i])*cos(thta1);
           dF3z[i]=(dM3[i]/radius[i])*cos(thta2);              
    } 
    

Log in to reply
 

京ICP备15017992号-2