粘性流体力学大作业

1、南京航空航天大学粘性流体力学大作业微型机翼设计报告一、题目及要求 某小型无人机重40kg，设计飞行速度100m/s，飞行高度2000m。使用Foil.html等课件作工具，设计其机翼。（1）应使该机翼在2度攻角时可产生足够升力保持飞机匀速平飞；（2）且尽量使附面层（尤其是上翼面）的压力梯度（或速度分布）不产生分离、或分离区尽量小；（3）分析估算摩擦阻力，应尽量减小摩阻。二、设计过程1、使用Foil.html等课件，设计其机翼。（1） 在完成公制单位等辅助设置后，选择指定的飞行速度，高度。（2） 在保持2度攻角情况下，设计机翼弯度、厚度，（3） 设计机翼弦长、翼展，（4） 利用输出功能分析机翼性

2、能及上下表面速度、压力等分布。2、结合机翼的表面压力（或速度）沿程分布，做2种以上方案进行对比分析，设计一个分离区尽量小的方案。3、利用Foil得到的机翼数据，分析估算摩擦阻力，应尽量减小摩阻。（1） 利用Foil得到的机翼数据，建立数据文件;（2） 编写附面层Karman积分计算的程序，读入你所设计机翼的数据，进行上下表面动量损失厚度的计算;（3） 附面层Karman积分计算采用以下湍流计算方法：其中无量纲参数和l满足：采用Thwaites方法：则当地摩阻为：根据F-S方程解和实验数据，可认为l和H都仅是的单变量函数，故得：将用表示的H和当地摩阻带入上式得：4步Runge-Kutta法步长示

3、意图解常微分方程的Runge-Kutta多步法：（4） 根据最后解得的附面层动量损失厚度计算机翼上下表面的摩擦阻力。（5） 利用整个计算分析系统，对不同设计方案的机翼开展摩擦阻力的对比分析。由计算得到的形状因子说明各个方案气流分离情况（以H>3.55为标准）。三、设计程序function OUTS=Drag_Airfoil% Generic time marching code solving the PDE for one dimensional wave :% Written by Huang Guoping, 2008/5/4nmax=19; % input the data of

4、 an airfoilDensity,Tem,Vupstream,Chord,Span,DataU,N_U,DataL,N_L=inputData();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.6'XL=Chord*DataL(:,1)' YL=Chord*DataL(:,2)' PL=DataL(:,3)*1000' VL=D

5、ataL(:,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:N_U dx(n) = dis(XU,YU,n); lengthU(n)= lengthU(n-1)+dx(n); if n=2 thetaU(n),HU(n)= BoundaryLayer_Flatplate(lengthU(n),

6、VU(n),Density,miu); else thetaU(n),HU(n)= BoundaryLayerEquation(dx(n),n,VU,Density,miu,thetaU(n-1); end %out=n, Density*VU(n)*lengthU(n)/miu/1e6, thetaU(n), HU(n)end % compute the boundary layer of airfoil's lower surfacelengthL(1)=0; thetaL(1)=0; CfL(1)=0; HL(1)=1;for n = 2:N_L dx(n) = dis(XU,Y

7、U,n); lengthL(n)= lengthL(n-1)+dx(n); if n=2 thetaL(n),HL(n)= BoundaryLayer_Flatplate(lengthL(n),VL(n),Density,miu); else thetaL(n),HL(n)= BoundaryLayerEquation(dx(n),n,VL,Density,miu,thetaL(n-1); end %out=n, Density*VL(n)*lengthL(n)/miu/1e6, thetaL(n), HL(n)end % plot the results of airfoil% output

8、 the Upper surface's parametersplotResults(lengthU,VU/Vsound,thetaU/(Chord*0.01),HU);% output the Lower surface's parametersplotResults(lengthL,VL/Vsound,thetaL/(Chord*0.01),HL); % Obtain the frictional DargDragU=thetaU(N_U)*Span*Density*Vupstream*Vupstream;DragL=thetaL(N_L)*Span*Density*Vup

9、stream*Vupstream;Drag =DragU+DragL; % END OF MAIN %function Density,Tem,Vupstream,Chord,Span,DataU,N_U,DataL,N_L=inputData()%N=input('enter no of grid points_');file1 = fopen('foil-0.dat', 'r');ccc=fscanf(file1, '%7f %7f %7f %7f %7f',5 1)'Density=ccc(1); Tem=ccc(2

10、); Vupstream=ccc(3); Chord=ccc(4); Span=ccc(5);tempc = fscanf(file1, '%25c',1 1);N_U = fscanf(file1, '%5i',1 1);tempc = fscanf(file1, '%25c',1 1);N_L = fscanf(file1, '%5i',1 1);fclose(file1); file2 = fopen('foil-U.dat', 'r');DataU = fscanf(file2, '

11、%8f %8f %7f %6f', 4 N_U)'fclose(file2); file3 = fopen('foil-L.dat', 'r');DataL = fscanf(file3, '%8f %8f %7f %6f', 4 N_L)'fclose(file3); % 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

12、)distance = sqrt(X(n)-X(n-1)2+(Y(n)-Y(n-1)2); % END function plotfoil(XU,YU,XL,YL)figurehold on;plot(XU,YU,'-o');plot(XL,YL,'-o');axis 'equal'hold off; % END function plotResults(L,V,theta,H)figureplot(L,V,'-d');figurehold on;plot(L,theta,'-d');plot(L,H,'-

13、o');hold off; % END function theta,H= BoundaryLayer_Flatplate(length,V,Density,miu)Rel =Density*V*length/miu; % Blasuis Solution for laminar flow theta =0.664*length/sqrt(Rel); Cf =0.664/sqrt(Rel); H =2.59; % % Algorithm for turbulent flow theta =0.0142*(Rel(6/7)*miu/(Density*V); Cf =0.026 *(Rel

14、(-1/7); H =1.375; % % END function LL,HH= Thwaites(Lamda) if Lamda>0.25 Lamda=0.25; else if Lamda<=-0.09; Lamda=-0.09; end end if Lamda>=0 LL=0.22+1.57*Lamda-1.8*Lamda2; HH=2.61-3.75*Lamda-5.24*Lamda2; else LL=0.22+1.042*Lamda+0.018*Lamda/(0.107+Lamda); HH=2.088+0.0731/(0.14+Lamda); end % E

15、ND function theta,HH= BoundaryLayerEquation(dx,n,V,Density,miu,theta1)% Solution of Runge-Kutta methodSita(1)=theta1;CF=0;HH=1;Sita(2)=Sita(1)+dx*(CF/2-(2+HH)*Sita(1)*(V(n)-V(n-1)/V(n)/dx)/8;for i=2:1:4 Lamda=Density*(Sita(i)2)*(V(n)-V(n-1)/miu/dx; LL,HH= Thwaites(Lamda); CF=2*miu*LL/Density/V(n)/Si

16、ta(i); Sita(i+1)=Sita(1)+dx*(CF/2-(2+HH)*Sita(i)*(V(n)-V(n-1)/V(n)/dx)/(2(4-i);end H=HH; Cf=CF; theta=Sita(i);% END function theta2,HH= BoundaryLayerEquation1(dx,n,V,Density,miu,theta1)% Solution of Euler method coupled with Iterationtheta2 = theta1; theta_error=theta1;Iter=0;while (theta_error>a

17、bs(theta1*0.001) && (Iter<10) Iter = Iter+1; theta_old= theta2; theta = 0.5*(theta1+theta2); ReTheta= Density*theta*0.5*(V(n)+V(n-1)/miu; Lamda = theta*ReTheta*(V(n)-V(n-1)/dx)/(0.5*(V(n)+V(n-1); Lamda = min(0.3,max(-0.09,Lamda); LL,HH= Thwaites(Lamda); FF=(LL-(2+HH)*Lamda)/ReTheta; theta

18、2= theta1 + dx*FF; theta2= max(theta1*0.01, theta2); theta_error =abs(theta2-theta_old);end % END 四、数据分析1、弯度不变，弦长和翼展变化：（1）第一组数据 :Camber = 3.0 % chord；Thickness = 9.0 % chord；Chord = 0.3 m；Span = 0.5 m。图中figure1-5分别为：a) 机翼外形；b) 上表面

19、的速度分布；c) 下表面的速度分布； d) 上表面H因子和Theta的变化；e) 下表面H因子和Theta的变化。详细输出数据： 上表面out =2.0000    9.1854    0.0000    0.0028    2.5900 out =3.0000    8.9624    0.0000

20、0;   0.0019    2.6659 out =4.0000    8.5760    0.0000    0.0013    2.7515 out =5.0000    8.2341    0.0000   &#

21、160;0.0011    2.7713 out = 6.0000    7.9220    0.0000    0.0009    2.8093out =7.0000    7.6396    0.0001    0.0007 

22、   2.8416out =8.0000    7.3572    0.0001    0.0006    2.9262 out =9.0000    7.0897    0.0001    0.0005    3.02

23、67 out =10.0000    6.8073    0.0001    0.0002    3.3481  out =11.0000    6.5397    0.0001    0.0001    3.4500 

24、out =12.0000    6.2722    0.0001    0.0001    3.5422  附面层出现分离 out =13.0000    6.0047    0.0001    0.0001    3.5918 附面层出现

25、分离 out =14.0000    5.7520    0.0002    0.0001    3.3212  附面层出现分离 out =15.0000    5.5142    0.0002    0.0001    

26、;3.5500 附面层出现分离out =16.0000    5.3061    0.0002    0.0001    3.5502  附面层出现分离 out =17.0000    5.1277    0.0002    0.0001  

27、0; 3.1988  附面层出现分离 out =18.0000    4.9643    0.0002    0.0001    3.7822  附面层出现分离 out =19.0000    3.8644    0.0005   

28、0;0.0000    3.6613 附面层出现分离下表面：out =2.0000    1.0553    0.0001    0.0082    2.5900 out =3.0000    2.5862    0.0000    0.0108&#

29、160;   1.9992 out =4.0000    4.0279    0.0000    0.0060    2.2220 out =5.0000    4.6670    0.0000    0.0035   

30、 2.4251 out =6.0000    4.9494    0.0000    0.0024    2.5082 out =7.0000    5.0683    0.0001    0.0018    2.5571

31、0;out =8.0000    5.0832    0.0001    0.0014    2.6019 out =9.0000    5.0386    0.0001    0.0012    2.6426 out =10.0

32、000    4.9791    0.0001    0.0010    2.6657 out =11.0000    4.9048    0.0001    0.0009    2.7014 out =12.0000  

33、;  4.8156    0.0001    0.0007    2.7571 out =13.0000    4.7413    0.0001    0.0007    2.7601 out =14.0000    4

34、.6967    0.0001    0.0007    2.7123 out =15.0000    4.6670    0.0001    0.0007    2.6907 out =16.0000    4.6670 

35、60;  0.0001    0.0007    2.6100 out =17.0000    4.6819    0.0001    0.0008    2.5503 out =18.0000    4.7413    

36、;0.0001    0.0011    3.6146 附面层出现分离 out =19.0000    3.8644    0.0002    0.0001    3.5523  附面层出现分离DragU =3.8960 DragL =1.8175 Drag 

37、=5.7134分析：以上数据表明，所设计机翼其上表面在第11个点就出现附面层分离，下表面只有最后两点出现附面层分离，因此所设计机翼需要进行改进。改进情况见后三组数据。（2）第二组数据：Camber = 3.0 % chord；Thickness = 9.0 % chord；Chord = 0.355 m；Span = 0.455 m。图中figure1-5分别为：a) 机翼外形；b) 上表面的速度分布；c) 下表面的速度分布； d) 上表面H因子和Theta的变化；e) 下表面H因子和Theta的变化。Drag= 5.