微型机翼设计报告
一、题目及要求
某小型无人机重40kg ，设计飞行速度100m/s，飞行高度2000m 。使用Foil.html 等课件作工具，设计其机翼。
（1）应使该机翼在2度攻角时可产生足够升力保持飞机匀速平飞；
（2）且尽量使附面层（尤其是上翼面）的压力梯度（或速度分布）不产生分离、或分离区尽量小；
（3）分析估算摩擦阻力，应尽量减小摩阻。
二、设计过程
1、使用Foil.html 等课件，设计其机翼。
（1） 在完成公制单位等辅助设置后，选择指定的飞行速度，高度。
（2） 在保持2度攻角情况下，设计机翼弯度、厚度，
（3） 设计机翼弦长、翼展，
（4） 利用输出功能分析机翼性能及上下表面速度、压力等分布。
2、结合机翼的表面压力（或速度）沿程分布，做2种以上方案进行对比分析，设计一个分离区尽量小的方案。
3、利用Foil 得到的机翼数据，分析估算摩擦阻力，应尽量减小摩阻。 （1） 利用Foil 得到的机翼数据，建立数据文件;
1 / 20
（2） 编写附面层Karman 积分计算的程序，读入你所设计机翼的数据，进行上下表面
动量损失厚度的计算;
d θC θdU e
=f -(2+H ) ⋅dx 2U e dx
（3） 附面层Karman 积分计算采用以下湍流计算方法：
其中无量纲参数λ和l 满足：
采用Thwaites 方法：
则当地摩阻为：
根据F-S 方程解和实验数据，可认为l 和H 都仅是λ的单变量函数，故得：
2 / 20
将用λ表示的H 和当地摩阻带入上式得：
解常微分方程的Runge-Kutta 多步法：
4步Runge-Kutta 法步长示意图
θ步1=θn +∆步1x ⋅R (θn , Cf n , H n ) θ步2=θn +∆步2x ⋅R (θ步1, Cf 步1, H 步1)
θ步m =θn +∆步m x ⋅R (θ步m-1, Cf 步m-1, H 步m-
1)
（4） 根据最后解得的附面层动量损失厚度θ计算机翼上下表面的摩擦阻力。
Drag 上=what (θ上, ρ∞, U ∞, ) ? Drag =Drag 上+Drag 下
（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/4 nmax=19;
% input the data of an airfoil
3 / 20
[Density,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 of airfoil plotfoil(XU,YU,XL,YL);
% compute the boundary layer of airfoil's upper surface lengthU(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)]=
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 surface lengthL(1)=0; thetaL(1)=0; CfL(1)=0; HL(1)=1; for n = 2:N_L
4 / 20
dx(n) = dis(XU,YU,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 the Upper surface's parameters
plotResults(lengthU,VU/Vsound,thetaU/(Chord*0.01),HU); % output the Lower surface's parameters
plotResults(lengthL,VL/Vsound,thetaL/(Chord*0.01),HL);
% Obtain the frictional Darg
DragU=thetaU(N_U)*Span*Density*Vupstream*Vupstream; DragL=thetaL(N_L)*Span*Density*Vupstream*Vupstream; Drag =DragU+DragL;
% END OF MAIN
5 / 20
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 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); 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, '%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);
6 / 20
% 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) figure hold on ;
plot(XU,YU,'-o' ); plot(XL,YL,'-o' ); axis 'equal' ; hold off ;
% END
function plotResults(L,V,theta,H) figure
plot(L,V,'-d' ); figure hold on ;
plot(L,theta,'-d' );
7 / 20
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))*miu/(Density*V); Cf =0.026 *(Rel^(-1/7)); H =1.375; %}
% END
function [LL,HH]= Thwaites(Lamda) if Lamda>0.25 Lamda=0.25;
8 / 20
if Lamda=0
LL=0.22+1.57*Lamda-1.8*Lamda^2; HH=2.61-3.75*Lamda-5.24*Lamda^2; else
LL=0.22+1.042*Lamda+0.018*Lamda/(0.107+Lamda); HH=2.088+0.0731/(0.14+Lamda); end
% END
function [theta,HH]=
BoundaryLayerEquation(dx,n,V,Density,miu,theta1) % Solution of Runge-Kutta method Sita(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)=Sita(1)+dx*(CF/2-(2+HH)*Sita(i)*(V(n)-V(n-1))/V(n)/dx)/(2^(4-i));
9 / 20
H=HH; Cf=CF;
theta=Sita(i); % END
function [theta2,HH]=
BoundaryLayerEquation1(dx,n,V,Density,miu,theta1) % Solution of Euler method coupled with Iteration theta2 = theta1; theta_error=theta1; Iter=0;
while (theta_error>abs(theta1*0.001)) && (Iter
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(theta1*0.01, theta2); theta_error =abs(theta2-theta_old); end
% END
10 / 20
四、数据分析
1、弯度不变，弦长和翼展变化：
（1）第一组数据:
Camber = 3.0 % chord；
Thickness = 9.0 % chord；
Chord = 0.3 m；
Span = 0.5 m。
11 / 20
详细输出数据：
上表面
out =2.0000 9.1854 0.0000 0.0028 2.5900
out =3.0000 8.9624 0.0000 0.0019 2.6659
out =4.0000 8.5760 0.0000 0.0013 2.7515
out =5.0000 8.2341 0.0000 0.0011 2.7713
out = 6.0000 7.9220 0.0000 0.0009 2.8093
12 / 20
out =8.0000 7.3572 0.0001 0.0006 2.9262
out =9.0000 7.0897 0.0001 0.0005 3.0267
out =10.0000 6.8073 0.0001 0.0002 3.3481
out =11.0000 6.5397 0.0001 0.0001 3.4500
out =12.0000 6.2722 0.0001 0.0001 3.5422—— 附面层出现分离
out =13.0000 6.0047 0.0001 0.0001 3.5918——附面层出现分离
out =14.0000 5.7520 0.0002 0.0001 3.3212 ——附面层出现分离
out =15.0000 5.5142 0.0002 0.0001 3.5500——附面层出现分离
out =16.0000 5.3061 0.0002 0.0001 3.5502 ——附面层出现分离
out =17.0000 5.1277 0.0002 0.0001 3.1988 ——附面层出现分离
out =18.0000 4.9643 0.0002 0.0001 3.7822 ——附面层出现分离
out =19.0000 3.8644 0.0005 0.0000 3.6613——附面层出现分离
下表面：
out =2.0000 1.0553 0.0001 0.0082 2.5900
out =3.0000 2.5862 0.0000 0.0108 1.9992
out =4.0000 4.0279 0.0000 0.0060 2.2220
13 / 20
out =6.0000 4.9494 0.0000 0.0024 2.5082
out =7.0000 5.0683 0.0001 0.0018 2.5571
out =8.0000 5.0832 0.0001 0.0014 2.6019
out =9.0000 5.0386 0.0001 0.0012 2.6426
out =10.0000 4.9791 0.0001 0.0010 2.6657
out =11.0000 4.9048 0.0001 0.0009 2.7014
out =12.0000 4.8156 0.0001 0.0007 2.7571
out =13.0000 4.7413 0.0001 0.0007 2.7601
out =14.0000 4.6967 0.0001 0.0007 2.7123
out =15.0000 4.6670 0.0001 0.0007 2.6907
out =16.0000 4.6670 0.0001 0.0007 2.6100
out =17.0000 4.6819 0.0001 0.0008 2.5503
out =18.0000 4.7413 0.0001 0.0011 3.6146 ——附面层出现分离
out =19.0000 3.8644 0.0002 0.0001 3.5523 ——附面层出现分离
DragU =3.8960
DragL =1.8175
14 / 20
Drag =5.7134
分析：以上数据表明，所设计机翼其上表面在第11个点就出现附面层分离，下表面只有最后两点出现附面层分离，因此所设计机翼需要进行改进。改进情况见后三组数据。
（2）第二组数据：
Camber = 3.0 % chord；
Thickness = 9.0 % chord；
Chord = 0.355 m；
Span = 0.455 m。
15 / 20
Drag= 5.3640
16 / 20
（3）结论：
比较上图可知，当camber 不变的时候，chord 和span 的变化对分离的点的影响较小，基本可以忽略不计。空气阻力却与chord 和span 有着密切的关系，当chord 增大，span 减小时，阻力会随之减小，反之则阻力随之增大。具体原因是由于机翼前缘附面层较薄，因此速度梯度较大，所以机翼前缘的粘性阻力较大，机翼沿流线方向向后则空气阻力随之减小。因此，机翼弦长较短，翼展较大时，相对的机翼前缘就比较长，所以空气阻力就较大，反之则空气阻力较小。
2、弯度变化，弦长和翼展变化忽略时：
（1）第一组数据：
Camber =1.0 % chord；
Thickness = 8.0 % chord；
Chord = 0.355 m；
Span = 0.7 m。
17 / 20
Drag =1.6908
（2）第二组数据：
Camber =5.0 % chord；
Thickness = 8.0 % chord；
Chord = 0.250 m；
Span = 0.500 m。
18 / 20
Drag =2.8527
19 / 20
（3）结论：
由上述过程中，可以发现当camber 增大，Thickness 相对不变时，其上表 面的分离位置向后移动，但是下表面的中部会出现分离点，因此结论是设计机翼时需要根据具体情况设计，并且需要进行多次的反复的修改和优化，以达到最优的设计。camber 增大，其摩擦阻力成增大趋势；根据第一二组数据可知：当chord 增大，span 减小时，其摩擦阻力成减小趋势。
五、心得体会
经过将近一星期的努力我终于独立完成了这次粘性流体力学的大作业！相信一开始大家都一样毫无头绪，但认真地看了一遍老师的大作业讲解PPT 后，才发现其实所有步骤老师都用图文精心地准备好了。知道了要做什么，我首先又恶补了一番
matlab ，把老师给的源程序一运行发现错误甚多。认真检查了才发现少了好多end ，函数都没结束. 。之后在老师留“？”的地方补好函数后运行就得到了5张图。在分析的过程中，我觉得自己对粘性流体力学这门课学地更加透彻了，同时也对老师留的
Foil.html 工具产生了浓厚的兴趣，甚至登陆了NASA 官网逛了一圈，收获颇丰。相信这次的大作业给大家都留下了深刻的印象!
20 / 20
微型机翼设计报告
一、题目及要求
某小型无人机重40kg ，设计飞行速度100m/s，飞行高度2000m 。使用Foil.html 等课件作工具，设计其机翼。
（1）应使该机翼在2度攻角时可产生足够升力保持飞机匀速平飞；
（2）且尽量使附面层（尤其是上翼面）的压力梯度（或速度分布）不产生分离、或分离区尽量小；
（3）分析估算摩擦阻力，应尽量减小摩阻。
二、设计过程
1、使用Foil.html 等课件，设计其机翼。
（1） 在完成公制单位等辅助设置后，选择指定的飞行速度，高度。
（2） 在保持2度攻角情况下，设计机翼弯度、厚度，
（3） 设计机翼弦长、翼展，
（4） 利用输出功能分析机翼性能及上下表面速度、压力等分布。
2、结合机翼的表面压力（或速度）沿程分布，做2种以上方案进行对比分析，设计一个分离区尽量小的方案。
3、利用Foil 得到的机翼数据，分析估算摩擦阻力，应尽量减小摩阻。 （1） 利用Foil 得到的机翼数据，建立数据文件;
1 / 20
（2） 编写附面层Karman 积分计算的程序，读入你所设计机翼的数据，进行上下表面
动量损失厚度的计算;
d θC θdU e
=f -(2+H ) ⋅dx 2U e dx
（3） 附面层Karman 积分计算采用以下湍流计算方法：
其中无量纲参数λ和l 满足：
采用Thwaites 方法：
则当地摩阻为：
根据F-S 方程解和实验数据，可认为l 和H 都仅是λ的单变量函数，故得：
2 / 20
将用λ表示的H 和当地摩阻带入上式得：
解常微分方程的Runge-Kutta 多步法：
4步Runge-Kutta 法步长示意图
θ步1=θn +∆步1x ⋅R (θn , Cf n , H n ) θ步2=θn +∆步2x ⋅R (θ步1, Cf 步1, H 步1)
θ步m =θn +∆步m x ⋅R (θ步m-1, Cf 步m-1, H 步m-
1)
（4） 根据最后解得的附面层动量损失厚度θ计算机翼上下表面的摩擦阻力。
Drag 上=what (θ上, ρ∞, U ∞, ) ? Drag =Drag 上+Drag 下
（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/4 nmax=19;
% input the data of an airfoil
3 / 20
[Density,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 of airfoil plotfoil(XU,YU,XL,YL);
% compute the boundary layer of airfoil's upper surface lengthU(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)]=
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 surface lengthL(1)=0; thetaL(1)=0; CfL(1)=0; HL(1)=1; for n = 2:N_L
4 / 20
dx(n) = dis(XU,YU,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 the Upper surface's parameters
plotResults(lengthU,VU/Vsound,thetaU/(Chord*0.01),HU); % output the Lower surface's parameters
plotResults(lengthL,VL/Vsound,thetaL/(Chord*0.01),HL);
% Obtain the frictional Darg
DragU=thetaU(N_U)*Span*Density*Vupstream*Vupstream; DragL=thetaL(N_L)*Span*Density*Vupstream*Vupstream; Drag =DragU+DragL;
% END OF MAIN
5 / 20
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 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); 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, '%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);
6 / 20
% 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) figure hold on ;
plot(XU,YU,'-o' ); plot(XL,YL,'-o' ); axis 'equal' ; hold off ;
% END
function plotResults(L,V,theta,H) figure
plot(L,V,'-d' ); figure hold on ;
plot(L,theta,'-d' );
7 / 20
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))*miu/(Density*V); Cf =0.026 *(Rel^(-1/7)); H =1.375; %}
% END
function [LL,HH]= Thwaites(Lamda) if Lamda>0.25 Lamda=0.25;
8 / 20
if Lamda=0
LL=0.22+1.57*Lamda-1.8*Lamda^2; HH=2.61-3.75*Lamda-5.24*Lamda^2; else
LL=0.22+1.042*Lamda+0.018*Lamda/(0.107+Lamda); HH=2.088+0.0731/(0.14+Lamda); end
% END
function [theta,HH]=
BoundaryLayerEquation(dx,n,V,Density,miu,theta1) % Solution of Runge-Kutta method Sita(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)=Sita(1)+dx*(CF/2-(2+HH)*Sita(i)*(V(n)-V(n-1))/V(n)/dx)/(2^(4-i));
9 / 20
H=HH; Cf=CF;
theta=Sita(i); % END
function [theta2,HH]=
BoundaryLayerEquation1(dx,n,V,Density,miu,theta1) % Solution of Euler method coupled with Iteration theta2 = theta1; theta_error=theta1; Iter=0;
while (theta_error>abs(theta1*0.001)) && (Iter
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(theta1*0.01, theta2); theta_error =abs(theta2-theta_old); end
% END
10 / 20
四、数据分析
1、弯度不变，弦长和翼展变化：
（1）第一组数据:
Camber = 3.0 % chord；
Thickness = 9.0 % chord；
Chord = 0.3 m；
Span = 0.5 m。
11 / 20
详细输出数据：
上表面
out =2.0000 9.1854 0.0000 0.0028 2.5900
out =3.0000 8.9624 0.0000 0.0019 2.6659
out =4.0000 8.5760 0.0000 0.0013 2.7515
out =5.0000 8.2341 0.0000 0.0011 2.7713
out = 6.0000 7.9220 0.0000 0.0009 2.8093
12 / 20
out =8.0000 7.3572 0.0001 0.0006 2.9262
out =9.0000 7.0897 0.0001 0.0005 3.0267
out =10.0000 6.8073 0.0001 0.0002 3.3481
out =11.0000 6.5397 0.0001 0.0001 3.4500
out =12.0000 6.2722 0.0001 0.0001 3.5422—— 附面层出现分离
out =13.0000 6.0047 0.0001 0.0001 3.5918——附面层出现分离
out =14.0000 5.7520 0.0002 0.0001 3.3212 ——附面层出现分离
out =15.0000 5.5142 0.0002 0.0001 3.5500——附面层出现分离
out =16.0000 5.3061 0.0002 0.0001 3.5502 ——附面层出现分离
out =17.0000 5.1277 0.0002 0.0001 3.1988 ——附面层出现分离
out =18.0000 4.9643 0.0002 0.0001 3.7822 ——附面层出现分离
out =19.0000 3.8644 0.0005 0.0000 3.6613——附面层出现分离
下表面：
out =2.0000 1.0553 0.0001 0.0082 2.5900
out =3.0000 2.5862 0.0000 0.0108 1.9992
out =4.0000 4.0279 0.0000 0.0060 2.2220
13 / 20
out =6.0000 4.9494 0.0000 0.0024 2.5082
out =7.0000 5.0683 0.0001 0.0018 2.5571
out =8.0000 5.0832 0.0001 0.0014 2.6019
out =9.0000 5.0386 0.0001 0.0012 2.6426
out =10.0000 4.9791 0.0001 0.0010 2.6657
out =11.0000 4.9048 0.0001 0.0009 2.7014
out =12.0000 4.8156 0.0001 0.0007 2.7571
out =13.0000 4.7413 0.0001 0.0007 2.7601
out =14.0000 4.6967 0.0001 0.0007 2.7123
out =15.0000 4.6670 0.0001 0.0007 2.6907
out =16.0000 4.6670 0.0001 0.0007 2.6100
out =17.0000 4.6819 0.0001 0.0008 2.5503
out =18.0000 4.7413 0.0001 0.0011 3.6146 ——附面层出现分离
out =19.0000 3.8644 0.0002 0.0001 3.5523 ——附面层出现分离
DragU =3.8960
DragL =1.8175
14 / 20
Drag =5.7134
分析：以上数据表明，所设计机翼其上表面在第11个点就出现附面层分离，下表面只有最后两点出现附面层分离，因此所设计机翼需要进行改进。改进情况见后三组数据。
（2）第二组数据：
Camber = 3.0 % chord；
Thickness = 9.0 % chord；
Chord = 0.355 m；
Span = 0.455 m。
15 / 20
Drag= 5.3640
16 / 20
（3）结论：
比较上图可知，当camber 不变的时候，chord 和span 的变化对分离的点的影响较小，基本可以忽略不计。空气阻力却与chord 和span 有着密切的关系，当chord 增大，span 减小时，阻力会随之减小，反之则阻力随之增大。具体原因是由于机翼前缘附面层较薄，因此速度梯度较大，所以机翼前缘的粘性阻力较大，机翼沿流线方向向后则空气阻力随之减小。因此，机翼弦长较短，翼展较大时，相对的机翼前缘就比较长，所以空气阻力就较大，反之则空气阻力较小。
2、弯度变化，弦长和翼展变化忽略时：
（1）第一组数据：
Camber =1.0 % chord；
Thickness = 8.0 % chord；
Chord = 0.355 m；
Span = 0.7 m。
17 / 20
Drag =1.6908
（2）第二组数据：
Camber =5.0 % chord；
Thickness = 8.0 % chord；
Chord = 0.250 m；
Span = 0.500 m。
18 / 20
Drag =2.8527
19 / 20
（3）结论：
由上述过程中，可以发现当camber 增大，Thickness 相对不变时，其上表 面的分离位置向后移动，但是下表面的中部会出现分离点，因此结论是设计机翼时需要根据具体情况设计，并且需要进行多次的反复的修改和优化，以达到最优的设计。camber 增大，其摩擦阻力成增大趋势；根据第一二组数据可知：当chord 增大，span 减小时，其摩擦阻力成减小趋势。
五、心得体会
经过将近一星期的努力我终于独立完成了这次粘性流体力学的大作业！相信一开始大家都一样毫无头绪，但认真地看了一遍老师的大作业讲解PPT 后，才发现其实所有步骤老师都用图文精心地准备好了。知道了要做什么，我首先又恶补了一番
matlab ，把老师给的源程序一运行发现错误甚多。认真检查了才发现少了好多end ，函数都没结束. 。之后在老师留“？”的地方补好函数后运行就得到了5张图。在分析的过程中，我觉得自己对粘性流体力学这门课学地更加透彻了，同时也对老师留的
Foil.html 工具产生了浓厚的兴趣，甚至登陆了NASA 官网逛了一圈，收获颇丰。相信这次的大作业给大家都留下了深刻的印象!
20 / 20