1、三.设计结果程序function OUTS=Drag_Airfoil% Example code for solving the Boundary 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=
2、Chord*DataU（:,2） PU=DataU（:,3）*1000 VU=DataU（:,4）/3.6XL=Chord*DataL（: YL=Chord*DataL（: PL=DataL（: VL=DataL（:% plot the shape of airfoilplotfoil（XU,YU,XL,YL）;% compute the boundary layer of airfoils upper surfacelengthU（1）=0; thetaU（1）=0; CfU（1）=0; HU（1）=1;for n = 2:nmax dx（n） = dis（XU,YU,n）; lengthU
3、（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）ends l
4、ower surfacelengthL（1）=0; thetaL（1）=0; CfL（1）=0; HL（1）=1; lengthL（n）= lengthL（n-1）+dx（n）; thetaL（n）,CfL（n）,HL（n）= BoundaryLayer_Flatplate（lengthL（n）,VL（n）,Density,miu）; thetaL（n）,CfL（n）,HL（n）= BoundaryLayerEquation（dx（n）,lengthL（n）,n,VL,Density,miu,thetaL（n-1）,CfL（n-1）,HL（n-1）; out=n, Density*VL（n）*
5、length（n）/miu/1e6, thetaL（n）, CfL（n）, HL（n）% 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/Vupstream,thetaU/（Chord*0.001）,CfU,HU）;plotResults（lengthL,VL/Vupstream,thetaL/（Chord*0.001）,CfL,HL）;DragU=
6、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=input（enter no of grid points_）;file1 = fopen（foil.dat, rccc=fscanf（file1, %7f %7f %7f %7f %7f,5 1）De
7、nsity=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, 4 nmax）DataL = fscanf（file1, fclose（file1）;% END function miu = Sutherland（Tem）miu0=1.4587e-6; Tem0=110.4;miu =miu0*（Tem）1.5）/（Tem+Tem0）;function distance=dis（
8、X,Y,n）distance = sqrt（X（n）-X（n-1）2+（Y（n）-Y（n-1）2）;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
9、 =0.0142*（Rel（6/7）*miu/（Density*V）; Cf =0.026 *（Rel（-1/7）; H =1.375;elsefunction Drag = DragP（nmax,X,Y,P）Drag = 0;for n=2: dy = Y（n）-Y（n-1）; Drag = Drag+0.5*dy*（P（n）+P（n-1）;function plotfoil（XU,YU,XL,YL）figurehold on;plot（XU,YU,-oplot（XL,YL,axis equalhold off;function theta,Cf,H= BoundaryLayerEquati
10、on（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.09; Lamuda（i）=-0.09;=0 l（i）=0.22+1.57*L
11、amuda（i）-1.8*Lamuda（i）2; Hh（i）=2.61-3.75*Lamuda（i）-5.24*Lamuda（i）2; 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）; CF（i）=2*miu*l（i）/Density/V（n）/Sita（i）; Sita（i+1）=Sita（1）+dx*（CF（i）/2-（2+Hh（i）*Sita（i）*（V（n）-V（n-1）/V（n）/dx）/（2（4-i）; if Lamuda（i）=-0.09
12、disp（附面层出现分离 H=Hh（i）; Cf=CF（i）; theta=Sita（i）;function plotResults（L,V,theta,Cf,H）plot（L,V,-rsplot（L,theta,%plot（L,Cf,-*plot（L,H,-+其中黑体部分是自编程序部分，具体是根据上文所述的内容而编出的。程序运行结果及结果数据分析见下文。四.设计结果分析1.关于分离点的位置影响因素：（1）、camber不变时，chord和span变化时 .第一组数据:Camber = 2.0 % chord , Thickness = 9.475 % chord , Chord = 0.35
13、9 m , Span = 0.868 m 机翼外形图像（1） （2） （3） （4）Drag=5.9134以下是输出详细数据：out =2.0000 9.1854 0.0000 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.8416ou
14、t =8.0000 7.3572 0.0001 0.0006 2.9262out =9.0000 7.0897 0.0001 0.0005 3.0267out =10.0000 6.8073 0.0001 0.0002 3.3481附面层出现分离out =11.0000 6.5397 0.0001 0.0001 3.5500out =12.0000 6.2722 0.0001 0.0001 3.5500out =13.0000 6.0047 0.0001 0.0001 3.5500out =14.0000 5.7520 0.0002 0.0001 3.5500out =15.0000 5.51
15、42 0.0002 0.0001 3.5500out =16.0000 5.3061 0.0002 0.0001 3.5500out =17.0000 5.1277 0.0002 0.0001 3.5500out =18.0000 4.9643 0.0002 0.0001 3.5500out =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.2220
16、out =5.0000 4.6670 0.0000 0.0035 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.000
17、1 0.0007 2.7571out =13.0000 4.7413 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.2146out =19.0000 3.8644 0.0002 0.0001 3.5500DragU
18、 =3.9960DragL =1.9175Drag =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）Drag=
19、 5.2640结论：由上图可知，当camber不变的时候，chord和span的变化对分离的点的影响较小，基本可以忽略不计。空气阻力却与chord和span有着密切的关系，当chord增大，span减小时，阻力会随之减小，反之则阻力随之增大。具体原因是由于机翼前缘附面层较薄，因此速度梯度较大，所以机翼前缘的粘性阻力较大，机翼沿流线方向向后则空气阻力随之减小。因此，机翼弦长较短，翼展较大时，相对的机翼前缘就比较长，所以空气阻力就较大，反之则空气阻力较小。（2）、camber变化时，chord和span的影响因素可忽略。.第三组数据：Camber = 0.4 % chord , Thickness
20、 = 6.5 % chord ,Chord = 0.365 m , Span = 0.982 m（3） （4）Drag =1.4508.第四组数据：Camber = 5.0 % chord , Thickness = 6.5 % chord ,Chord = 0.31 m , Span = 0.64 mDrag =2.1672.第五组数据：Camber = 9.5 % chord , Thickness = 6.5 % chord ,Chord = 0.256 m , Span = 0.525 mDrag = 3.7562由上述过程中，可以发现当camber增大，Thickness相对不变时，
21、其上表面的分离位置向后移动，但是下表面的中部会出现分离点，因此结论是设计机翼时需要根据具体情况设计，并且需要进行多次的反复的修改和优化，以达到最优的设计。camber增大，其摩擦阻力成增大趋势；根据第一二组数据可知：当chord增大，span减小时，其摩擦阻力成减小趋势。五.机翼设计心得经过长时间的努力和坚持，我终于完成了粘性流体力学的大作业。粘性流体力学大作业是让我们将理论的学习应用于实践，并且让它得到升华，让我对粘性流体力学这门课程有了更加深刻的了解，并且感受到了这门课程的独特和无穷的魅力。本次的大作业是锻炼我们实际操作的能力，培养了我们工程设计的理念。首先，在刚开始设计时，我只是在Foi
22、l.html这个软件下凭借自己的主观感觉进行设计，后来通过Matlab程序运行计算出结果后我通过多组数据进行对比，采用了变量的方法，发现了设计机翼中的那种潜在的规律，这对于我学习和理解粘性流体力学这门课程有着很大的帮助。其次， 在实际设计过程中，我们不仅需要掌握粘性流体力学的基本知识，还要熟练运用MATLAB进行数据分析，因此Matlab编程就是很关键的一个环节。由于原来接触的Matlab只是简单的几句语句，因此，这次编程对我来说是有着很大的困难，但是经过不懈的努力和与同学们的讨论，我终于编出了用于解决这个问题的程序来，这大大加深了我对龙格库塔差值方法的理解和掌握，以及认识了它的应用。在进行大作业的整个过程中需要各方面的知识的融会贯通，才能进行初步的程序设计。然后，在数据分析的阶段，我设计了很多组机翼，分析每个参数对机翼的影响，因此我耗费了大量的时间，但是我从中也学到了许多东西，加深了在课堂上学习内容的理解。这次机翼的设计对我来说可以认为是一种挑战，但它真正的培养了我们所需的工程设计的基本本领，让我们学到了很多课堂上无法学到了知识。
