粘性流體力學(xué)大作業(yè)

1、南京航空航天大學(xué)粘性流體力學(xué)大作業(yè)微型機(jī)翼設(shè)計(jì)報(bào)告一、題目及要求 某小型無(wú)人機(jī)重40kg，設(shè)計(jì)飛行速度100m/s，飛行高度2000m。使用foil.html等課件作工具，設(shè)計(jì)其機(jī)翼。（1）應(yīng)使該機(jī)翼在2度攻角時(shí)可產(chǎn)生足夠升力保持飛機(jī)勻速平飛；（2）且盡量使附面層（尤其是上翼面）的壓力梯度（或速度分布）不產(chǎn)生分離、或分離區(qū)盡量??；（3）分析估算摩擦阻力，應(yīng)盡量減小摩阻。二、設(shè)計(jì)過(guò)程1、使用foil.html等課件，設(shè)計(jì)其機(jī)翼。（1） 在完成公制單位等輔助設(shè)置后，選擇指定的飛行速度，高度。（2） 在保持2度攻角情況下，設(shè)計(jì)機(jī)翼彎度、厚度，（3） 設(shè)計(jì)機(jī)翼弦長(zhǎng)、翼展，（4） 利用輸出功能分析機(jī)翼性

2、能及上下表面速度、壓力等分布。2、結(jié)合機(jī)翼的表面壓力（或速度）沿程分布，做2種以上方案進(jìn)行對(duì)比分析，設(shè)計(jì)一個(gè)分離區(qū)盡量小的方案。3、利用foil得到的機(jī)翼數(shù)據(jù)，分析估算摩擦阻力，應(yīng)盡量減小摩阻。（1） 利用foil得到的機(jī)翼數(shù)據(jù)，建立數(shù)據(jù)文件;（2） 編寫附面層karman積分計(jì)算的程序，讀入你所設(shè)計(jì)機(jī)翼的數(shù)據(jù)，進(jìn)行上下表面動(dòng)量損失厚度的計(jì)算;（3） 附面層karman積分計(jì)算采用以下湍流計(jì)算方法：其中無(wú)量綱參數(shù)和l滿足：采用thwaites方法：則當(dāng)?shù)啬ψ铻椋焊鶕?jù)f-s方程解和實(shí)驗(yàn)數(shù)據(jù)，可認(rèn)為l和h都僅是的單變量函數(shù)，故得：將用表示的h和當(dāng)?shù)啬ψ鑾肷鲜降茫?步runge-kutta法步長(zhǎng)示

3、意圖解常微分方程的runge-kutta多步法：（4） 根據(jù)最后解得的附面層動(dòng)量損失厚度計(jì)算機(jī)翼上下表面的摩擦阻力。（5） 利用整個(gè)計(jì)算分析系統(tǒng)，對(duì)不同設(shè)計(jì)方案的機(jī)翼開展摩擦阻力的對(duì)比分析。由計(jì)算得到的形狀因子說(shuō)明各個(gè)方案氣流分離情況（以h3.55為標(biāo)準(zhǔn)）。三、設(shè)計(jì)程序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 四、數(shù)據(jù)分析1、彎度不變，弦長(zhǎng)和翼展變化：（1）第一組數(shù)據(jù) :camber=3.0%chord；thickness=9.0%chord；chord=0.3m；span=0.5m。圖中figure1-5分別為：a) 機(jī)翼外形；b) 上表面的速度分布；c) 下表面的速度分布；d) 上表面h因子和theta的變化；e) 下表面h因子和theta的變化。詳細(xì)輸出數(shù)據(jù)：上表面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 附面層出現(xiàn)分離out=13.00006.00470.00010.00013.5918 附面層出現(xiàn)分離out=14.00005.75200.00020.00013.3212附面層出現(xiàn)分離out=15.00005.51420.00020.00013.5500 附面層出現(xiàn)分離out=16.00005.30610.00020.00013.5502附面層出現(xiàn)分離out=17.00005.12770.00020.00013.1988附面層出現(xiàn)分離out=18.00004.96430.00020.00013.7822附面層出現(xiàn)分離out=19.00003

20、.86440.00050.00003.6613 附面層出現(xiàn)分離下表面：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附面層出現(xiàn)分離out=19.00003.86440.00020.00013.5523附面層出現(xiàn)分離dragu=3.8960dragl=1.8175drag=5.7134分析：以上數(shù)據(jù)表明，所設(shè)計(jì)機(jī)翼其上表面在第11個(gè)點(diǎn)就出現(xiàn)附面層分離，下表面只有最后兩點(diǎn)出現(xiàn)附面層分離，因此所設(shè)計(jì)機(jī)翼需要進(jìn)行改進(jìn)。改進(jìn)情況見后三組數(shù)據(jù)。（2）第二組數(shù)據(jù)：camber=3.0%chord；thickness=9.0%chord；chord=0.355m；span=0.455m。圖中figure1-5分別為：a) 機(jī)翼外形；b) 上表面的速度分布；c) 下表面的速度分布；

23、d) 上表面h因子和theta的變化；e) 下表面h因子和theta的變化。drag=5.3640（3）結(jié)論：比較上圖可知，當(dāng)camber不變的時(shí)候，chord和span的變化對(duì)分離的點(diǎn)的影響較小，基本可以忽略不計(jì)?？諝庾枇s與chord和span有著密切的關(guān)系，當(dāng)chord增大，span減小時(shí)，阻力會(huì)隨之減小，反之則阻力隨之增大。具體原因是由于機(jī)翼前緣附面層較薄，因此速度梯度較大，所以機(jī)翼前緣的粘性阻力較大，機(jī)翼沿流線方向向后則空氣阻力隨之減小。因此，機(jī)翼弦長(zhǎng)較短，翼展較大時(shí)，相對(duì)的機(jī)翼前緣就比較長(zhǎng)，所以空氣阻力就較大，反之則空氣阻力較小。2、彎度變化，弦長(zhǎng)和翼展變化忽略時(shí)：（1）第一組數(shù)據(jù)：camber=1.0%chord；thickness=8.0%