Posts made by maguolin

部分程序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);              
} 

部分程序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;
 } 

部分程序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);
    }                
} 

各位老师好，我最近按照中科院王胜军博士的博士论文编写了一个风力机致动线模型的UDF，但一计算要么发散，要么浮点溢出。
我的致动线的程序，我是把叶片分成40个叶素，先确每个叶素的中心点，风轮初始位置及运行的每一时刻，叶素中心的坐标好确定，但这个坐标不一定的网格中心的坐标，此时需要找到与之最近的网格作为线上的网格，并确定其中心坐标及相应网格中心的速度，然后返回去用BEM求解相应网格点的体积力，再用高斯分布分散到此网格周围，现在的问题的与叶素中心最近的那个网格怎么找呢？我用的一个遍历，可能不太对，还有我那个体积力源项可能也不太对，并行计算时，一开始就浮点溢出，串行时一开始就发散，我原来以为可能计算量大浮点溢出，但为什么串行也不行呢？所以我后来想着是我的UDDF出了问题，但是，我实在找不出问题出在哪里了。
可能这个地方处理的不好

麻烦各位老师帮我看看，有没有更好一些的程序实现。
程序流程图：

这是部分程序代码，全部的好像传不上去，比较多。

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)   
     }
  }