粘性流体力学大作业

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%