粘性流体力学大作业

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） 利用整个计算分析系统，对不同设计方案的机翼开展摩擦阻力的对比分析。由计算得到的形状因子说明各个方案气流分离情况（以H3.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 an

4、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=DataL(:,4)/3.6; % plot the shape

5、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: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),VU(n),Density,miu); else thetaU(n),HU(n)=

6、 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 airfoils lower surfacelengthL(1)=0; thetaL(1)=0; CfL(1)=0; HL(1)=1;for n = 2:N_L dx(n) = dis(XU,YU,n); lengthL(n)= lengthL(n-1)+dx(n); if n=2 t

7、hetaL(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 the Upper surfaces parametersplotResults(leng

8、thU,VU/Vsound,thetaU/(Chord*0.01),HU);% output the Lower surfaces 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*Vupstream*Vupstream;Drag =DragU+DragL; % END OF MAIN %funct

9、ion 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); Vupstream=ccc(3); Chord=ccc(4); Span=ccc(5);tempc = fscanf(file1, %25c,1 1);N_U = fscanf(file1, %

10、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, %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)

11、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)-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);fig

12、urehold on;plot(L,theta,-d);plot(L,H,-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)

13、*miu/(Density*V); Cf =0.026 *(Rel(-1/7); H =1.375; % % END function LL,HH= Thwaites(Lamda) if Lamda0.25 Lamda=0.25; else 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 % END function thet

14、a,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)/Sita(i); Sita(i+1)

15、=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_errorabs(theta1*0.001) & (

16、Iter10) 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; theta2= theta1 + dx*FF; theta2= max(th

17、eta1*0.01, theta2); theta_error =abs(theta2-theta_old);end % END 四、数据分析1、弯度不变，弦长和翼展变化：（1）第一组数据 :Camber=3.0%chord；Thickness=9.0%chord；Chord=0.3m；Span=0.5m。图中figure1-5分别为：a) 机翼外形；b) 上表面的速度分布；c) 下表面的速度分布；d) 上表面H因子和Theta的变化；e) 下表面H因子和Theta的变化。详细输出数据：上表面out=2.00009.18540.00000.00282.5900out=3.00008.96240

18、.00000.00192.6659out=4.00008.57600.00000.00132.7515out=5.00008.23410.00000.00112.7713out=6.00007.92200.00000.00092.8093out=7.00007.63960.00010.00072.8416out=8.00007.35720.00010.00062.9262out=9.00007.08970.00010.00053.0267out=10.00006.80730.00010.00023.3481out=11.00006.53970.00010.00013.4500out=12.00

19、006.27220.00010.00013.5422 附面层出现分离out=13.00006.00470.00010.00013.5918 附面层出现分离out=14.00005.75200.00020.00013.3212附面层出现分离out=15.00005.51420.00020.00013.5500 附面层出现分离out=16.00005.30610.00020.00013.5502附面层出现分离out=17.00005.12770.00020.00013.1988附面层出现分离out=18.00004.96430.00020.00013.7822附面层出现分离out=19.00003

20、.86440.00050.00003.6613 附面层出现分离下表面：out=2.00001.05530.00010.00822.5900out=3.00002.58620.00000.01081.9992out=4.00004.02790.00000.00602.2220out=5.00004.66700.00000.00352.4251out=6.00004.94940.00000.00242.5082out=7.00005.06830.00010.00182.5571out=8.00005.08320.00010.00142.6019out=9.00005.03860.00010.001

21、22.6426out=10.00004.97910.00010.00102.6657out=11.00004.90480.00010.00092.7014out=12.00004.81560.00010.00072.7571out=13.00004.74130.00010.00072.7601out=14.00004.69670.00010.00072.7123out=15.00004.66700.00010.00072.6907out=16.00004.66700.00010.00072.6100out=17.00004.68190.00010.00082.5503out=18.00004.

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

23、d) 上表面H因子和Theta的变化；e) 下表面H因子和Theta的变化。Drag=5.3640（3）结论：比较上图可知，当camber不变的时候，chord和span的变化对分离的点的影响较小，基本可以忽略不计。空气阻力却与chord和span有着密切的关系，当chord增大，span减小时，阻力会随之减小，反之则阻力随之增大。具体原因是由于机翼前缘附面层较薄，因此速度梯度较大，所以机翼前缘的粘性阻力较大，机翼沿流线方向向后则空气阻力随之减小。因此，机翼弦长较短，翼展较大时，相对的机翼前缘就比较长，所以空气阻力就较大，反之则空气阻力较小。2、弯度变化，弦长和翼展变化忽略时：（1）第一组数据：Camber=1.0%chord；Thickness=8.0%