微型飞机机翼设计大作业

1、粘性流体力学大作业报告能源与动力学院微型飞机机翼设计报告学号： 姓名： 2010-4-27微型飞机机翼设计报告一.设计题目及要求某小型无人机重100kg，设计飞行速度100m/s，飞行高度3000m。使用Foil.html等课件作工具，设计其机翼，（1）应使该机翼在5度攻角时可产生足够升力保持飞机匀速平飞，（2）且尽量使附面层（尤其是上翼面）的压力梯度（或速度分布）不产生分离、或分离区尽量小；（3）分析估算摩擦阻力，应尽量减小摩阻。二.设计过程（1）使用Foil.html等课件，设计其机翼。（2）利用Foil得到的机翼数据，分析估算摩擦阻力，应尽量减小摩阻。1、利用Foil得到的机翼数据，建立

2、数据文件；2、编写附面层Karman积分计算的程序，读入你所设计机翼的数据，进行上下表面动量损失厚度的计算源公式：其中为动量损失厚度.在此MATLAB程序中我们采用简化算法即 其中：采用Thwaites方法：计算的取值范围为-0.09，+0.25 若计算中出现> +0.25，则取为+0.25若计算中出现< - 0.09，则取为- 0.09为了计算此积分，我们采用了龙格库塔的积分方法，其积分方法如下所示 其中：三.设计结果程序150207105班 粘性流体力学大作业报告function OUTS=Drag_Airfoil% Example code for solving the B

3、oundary Layer of airfoil% Written by Huang Guoping, 2007/4/10nmax=19; % input the data of an airfoilDensity,Tem,Vupstream,Chord,Span,DataU,DataL=inputData(nmax);miu = Sutherland(Tem); Vsound=sqrt(1.4*287.2*Tem);XU=Chord*DataU(:,1)' YU=Chord*DataU(:,2)' PU=DataU(:,3)*1000' VU=DataU(:,4)/3

4、.6'XL=Chord*DataL(:,1)' YL=Chord*DataL(:,2)' PL=DataL(:,3)*1000' VL=DataL(:,4)/3.6'% plot the shape of airfoilplotfoil(XU,YU,XL,YL);% compute the boundary layer of airfoil's upper surfacelengthU(1)=0; thetaU(1)=0; CfU(1)=0; HU(1)=1;for n = 2:nmax dx(n) = dis(XU,YU,n); lengthU

5、(n)= lengthU(n-1)+dx(n); if n=2 thetaU(n),CfU(n),HU(n)= BoundaryLayer_Flatplate(lengthU(n),VU(n),Density,miu); else thetaU(n),CfU(n),HU(n)= BoundaryLayerEquation(dx(n),lengthU(n),n,VU,Density,miu,thetaU(n-1),CfU(n-1),HU(n-1); end out=n, Density*VU(n)*length(n)/miu/1e6, thetaU(n), CfU(n), HU(n)end% c

6、ompute the boundary layer of airfoil's lower surfacelengthL(1)=0; thetaL(1)=0; CfL(1)=0; HL(1)=1;for n = 2:nmax dx(n) = dis(XU,YU,n); lengthL(n)= lengthL(n-1)+dx(n); if n=2 thetaL(n),CfL(n),HL(n)= BoundaryLayer_Flatplate(lengthL(n),VL(n),Density,miu); else thetaL(n),CfL(n),HL(n)= BoundaryLayerEq

7、uation(dx(n),lengthL(n),n,VL,Density,miu,thetaL(n-1),CfL(n-1),HL(n-1); end out=n, Density*VL(n)*length(n)/miu/1e6, thetaL(n), CfL(n), HL(n)end% compute the Pressure drag DragPU = DragP(nmax,XU,YU,PU)*Span;DragPL =-DragP(nmax,XL,YL,PL)*Span;% plot the results of airfoilplotResults(lengthU,VU/Vupstrea

8、m,thetaU/(Chord*0.001),CfU,HU);plotResults(lengthL,VL/Vupstream,thetaL/(Chord*0.001),CfL,HL);DragU=thetaU(nmax)*Span*Density*Vupstream*VupstreamDragL=thetaL(nmax)*Span*Density*Vupstream*VupstreamDrag =DragU+DragL% END OF MAIN%function Density,Tem,Vupstream,Chord,Span,DataU,DataL=inputData(nmax)%N=in

9、put('enter no of grid points_');file1 = fopen('foil.dat', 'r');ccc=fscanf(file1, '%7f %7f %7f %7f %7f',5 1)'Density=ccc(1); Tem=ccc(2); Vupstream=ccc(3); Chord=ccc(4); Span=ccc(5);tempc = fscanf(file1, '%20c',1 1);DataU = fscanf(file1, '%8f %8f %7f %6f

10、', 4 nmax)'tempc = fscanf(file1, '%20c',1 1);DataL = fscanf(file1, '%8f %8f %7f %6f', 4 nmax)'fclose(file1);% END function miu = Sutherland(Tem)miu0=1.4587e-6; Tem0=110.4;miu =miu0*(Tem)1.5)/(Tem+Tem0); % END function distance=dis(X,Y,n)distance = sqrt(X(n)-X(n-1)2+(Y(n)-

11、Y(n-1)2);% END function theta,Cf,H= BoundaryLayer_Flatplate(length,V,Density,miu)Rel =Density*V*length/miu; II =1;if II=1 % Blasuis Solution for laminar flow theta =0.664*length/sqrt(Rel); Cf =0.664/sqrt(Rel); H =2.59;elseif II=2 % Algorithm for turbulent flow theta =0.0142*(Rel(6/7)*miu/(Density*V)

12、; Cf =0.026 *(Rel(-1/7); H =1.375;elseend% END function Drag = DragP(nmax,X,Y,P)Drag = 0;for n=2:nmax dy = Y(n)-Y(n-1); Drag = Drag+0.5*dy*(P(n)+P(n-1); end% END function plotfoil(XU,YU,XL,YL)figurehold on;plot(XU,YU,'-o');plot(XL,YL,'-o');axis 'equal'hold off;% END function

13、theta,Cf,H= BoundaryLayerEquation(dx,length,n,V,Density,miu,theta1,Cf1,H1)Sita(1)=theta1;CF(1)=Cf1;Hh(1)=H1;Sita(2)=Sita(1)+dx*(CF(1)/2-(2+Hh(1)*Sita(1)*(V(n)-V(n-1)/V(n)/dx)/8;for i=2:1:4 Lamuda(i)=Density*(Sita(i)2)*(V(n)-V(n-1)/miu/dx; if Lamuda(i)>0.25 Lamuda(i)=0.25; elseif Lamuda(i)<=-0.

14、09; Lamuda(i)=-0.09; end if Lamuda(i)>=0 l(i)=0.22+1.57*Lamuda(i)-1.8*Lamuda(i)2; Hh(i)=2.61-3.75*Lamuda(i)-5.24*Lamuda(i)2; else l(i)=0.22+1.042*Lamuda(i)+0.018*Lamuda(i)/(0.107+Lamuda(i); Hh(i)=2.088+0.0731/(0.14+Lamuda(i); end CF(i)=2*miu*l(i)/Density/V(n)/Sita(i); Sita(i+1)=Sita(1)+dx*(CF(i)/

15、2-(2+Hh(i)*Sita(i)*(V(n)-V(n-1)/V(n)/dx)/(2(4-i);end if Lamuda(i)=-0.09 disp('附面层出现分离'); end H=Hh(i); Cf=CF(i); theta=Sita(i);% END function plotResults(L,V,theta,Cf,H)figurehold on;plot(L,V,'-rs');hold off;figurehold on;plot(L,theta,'-o');hold off;%plot(L,Cf,'-*');ho

16、ld on;plot(L,H,'-+');hold off;% END 其中黑体部分是自编程序部分，具体是根据上文所述的内容而编出的。程序运行结果及结果数据分析见下文。四.设计结果分析1.关于分离点的位置影响因素：（1）、camber不变时，chord和span变化时 .第一组数据:Camber = 2.0 % chord , Thickness = 9.475 % chord , Chord = 0.359 m , Span = 0.868 m 机翼外形图像(1) (2) (3) (4)Drag=5.9134以下是输出详细数据：out =2.0000 9.1854 0.000

17、0 0.0028 2.5900out =3.0000 8.9624 0.0000 0.0019 2.6659out =4.0000 8.5760 0.0000 0.0013 2.7515out =5.0000 8.2341 0.0000 0.0011 2.7713out = 6.0000 7.9220 0.0000 0.0009 2.8093out =7.0000 7.6396 0.0001 0.0007 2.8416out =8.0000 7.3572 0.0001 0.0006 2.9262out =9.0000 7.0897 0.0001 0.0005 3.0267out =10.000

18、0 6.8073 0.0001 0.0002 3.3481附面层出现分离out =11.0000 6.5397 0.0001 0.0001 3.5500附面层出现分离out =12.0000 6.2722 0.0001 0.0001 3.5500附面层出现分离out =13.0000 6.0047 0.0001 0.0001 3.5500附面层出现分离out =14.0000 5.7520 0.0002 0.0001 3.5500附面层出现分离out =15.0000 5.5142 0.0002 0.0001 3.5500附面层出现分离out =16.0000 5.3061 0.0002 0.

19、0001 3.5500附面层出现分离out =17.0000 5.1277 0.0002 0.0001 3.5500附面层出现分离out =18.0000 4.9643 0.0002 0.0001 3.5500附面层出现分离out =19.0000 3.8644 0.0005 0.0000 3.5500out =2.0000 1.0553 0.0001 0.0082 2.5900out =3.0000 2.5862 0.0000 0.0108 1.9992out =4.0000 4.0279 0.0000 0.0060 2.2220out =5.0000 4.6670 0.0000 0.003

20、5 2.4251out =6.0000 4.9494 0.0000 0.0024 2.5082out =7.0000 5.0683 0.0001 0.0018 2.5571out =8.0000 5.0832 0.0001 0.0014 2.6019out =9.0000 5.0386 0.0001 0.0012 2.6426out =10.0000 4.9791 0.0001 0.0010 2.6657out =11.0000 4.9048 0.0001 0.0009 2.7014out =12.0000 4.8156 0.0001 0.0007 2.7571out =13.0000 4.7

21、413 0.0001 0.0007 2.7601out =14.0000 4.6967 0.0001 0.0007 2.7123out =15.0000 4.6670 0.0001 0.0007 2.6907out =16.0000 4.6670 0.0001 0.0007 2.6100out =17.0000 4.6819 0.0001 0.0008 2.5503out =18.0000 4.7413 0.0001 0.0011 2.2146附面层出现分离out =19.0000 3.8644 0.0002 0.0001 3.5500DragU =3.9960DragL =1.9175Dra

22、g =5.9134注：上图中 （1）代表上表面的速度分布； （2）代表下表面的速度分布； （3）代表上表面H因子和Theta的变化； （4）代表下表面H因子和Theta的变化。 下列图像规则同上。分析：以上数据表明，所设计机翼其上表面在第9个点出现附面层分离，下表面只有最后一点出现附面层分离，因此所设计机翼需要进行改进。具体改进参见数据3、4、5.第二组数据:Camber = 2.0 % chord , Thickness = 9.475 % chord ,Chord = 0.377 m , Span = 0.754 m (1) (2) （3） （4）Drag= 5.2640结论：由上图可知，

23、当camber不变的时候，chord和span的变化对分离的点的影响较小，基本可以忽略不计。空气阻力却与chord和span有着密切的关系，当chord增大，span减小时，阻力会随之减小，反之则阻力随之增大。具体原因是由于机翼前缘附面层较薄，因此速度梯度较大，所以机翼前缘的粘性阻力较大，机翼沿流线方向向后则空气阻力随之减小。因此，机翼弦长较短，翼展较大时，相对的机翼前缘就比较长，所以空气阻力就较大，反之则空气阻力较小。（2）、camber变化时，chord和span的影响因素可忽略。.第三组数据：Camber = 0.4 % chord , Thickness = 6.5 % chord ,Chord = 0.365 m , Span = 0.982 m（1） （2）（3） （4）Drag =1.4508.第四组数据：Camber = 5.0 % chord , Thickness = 6.5 % chord ,Chord = 0.31 m , Span = 0.64 m（1） (2)（3） （4）Drag =2.1672.第五组数据：Camber = 9.5 % chord , Thickness = 6.5 % chord ,Chord = 0.256 m , Span = 0.525 m（1） （2）（3） （4）Drag = 3.7562结论：由上述过程中，可以发