Navigation

    CFD中文网

    CFD中文网

    • Login
    • Search
    • 最新

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

    Fluent
    3
    7
    1433
    Loading More Posts
    • Oldest to Newest
    • Newest to Oldest
    • Most Votes
    Reply
    • Reply as topic
    Log in to reply
    This topic has been deleted. Only users with topic management privileges can see it.
    • M
      maguolin last edited by

      各位老师好,我最近按照中科院王胜军博士的博士论文编写了一个风力机致动线模型的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)   
           }
        }
      
      1 Reply Last reply Reply Quote
      • M
        maguolin last edited by

        部分程序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);
            }                
        } 
        
        1 Reply Last reply Reply Quote
        • M
          maguolin last edited by

          部分程序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;
           } 
          
          1 Reply Last reply Reply Quote
          • M
            maguolin last edited by

            部分程序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);              
            } 
            
            五好青年 1 Reply Last reply Reply Quote
            • D
              Darry6682 last edited by

              你好请问还在研究AL+CFD吗? 可以交流一下吗

              五好青年 1 Reply Last reply Reply Quote
              • 五好青年
                五好青年 @Darry6682 last edited by

                @Darry6682 我一直做致动线的很多年,你有问题可以交流,我Email:jrw1992@163.com

                I am a CFD machine with no emotions. Welcome to browse my Bilibili, search: seeeeeeeeeeer

                1 Reply Last reply Reply Quote
                • 五好青年
                  五好青年 @maguolin last edited by

                  @maguolin
                  嗨,国林。能邮件交流下致动线吗?jrw1992@163.com
                  期待和你探讨

                  I am a CFD machine with no emotions. Welcome to browse my Bilibili, search: seeeeeeeeeeer

                  1 Reply Last reply Reply Quote
                  • First post
                    Last post

                  CFD中文网 | 东岳流体 | 京ICP备15017992号-2