关于用fluent UDF 编写风力机致动线模型
-
各位老师好，我最近按照中科院王胜军博士的博士论文编写了一个风力机致动线模型的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) } } -
部分程序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); } -
@Darry6682 好呀，我的邮箱是 maguolin1234@126.com
