粘性流体力学大作业.docx
《粘性流体力学大作业.docx》由会员分享，可在线阅读，更多相关《粘性流体力学大作业.docx（21页珍藏版）》请在冰点文库上搜索。
粘性流体力学大作业
微型机翼设计报告
一、题目及要求
某小型无人机重40kg，设计飞行速度100m/s，飞行高度2000m。
使用Foil.html等课件作工具，设计其机翼。
（1）应使该机翼在2度攻角时可产生足够升力保持飞机匀速平飞；
（2）且尽量使附面层（尤其是上翼面）的压力梯度（或速度分布）不产生分离、或分离区尽量小；
（3）分析估算摩擦阻力，应尽量减小摩阻。
1、设计过程
1、使用Foil.html等课件，设计其机翼。
（1）在完成公制单位等辅助设置后，选择指定的飞行速度，高度。
（2）在保持2度攻角情况下，设计机翼弯度、厚度，
（3）设计机翼弦长、翼展，
（4）利用输出功能分析机翼性能及上下表面速度、压力等分布。
2、结合机翼的表面压力（或速度）沿程分布，做2种以上方案进行对比分析，设计一个
分离区尽量小的方案。
3、利用Foil得到的机翼数据，分析估算摩擦阻力，应尽量减小摩阻。
（1）利用Foil得到的机翼数据，建立数据文件；
（2）编写附面层Karman积分计算的程序，读入你所设计机翼的数据，进行上下表面
动量损失厚度的计算
（3）附面层Karman积分计算采用以下湍流计算方法:
其中无量纲参数入和l满足:
l（）,H
H（）
采用Thwaites方法:
Ue
y2
2
2
y0:
udP-警ydxdx
dx
则当地摩阻为：
d1
dxRe
l
Re
0:
l0.22
1.57
1.82,
H
2.613.755.24
0.018
c—0.0731
0:
l0.22
1.402
H
2.088
0.107
0.14
0.09:
出现分离
入的单变量函数，故得:
根据F-S方程解和实验数据，可认为丨和H都仅是
0.09
将用入表示的H和当地摩阻带入上式得:
（2H）
Ue
（2H）—-
Re
1
Re
[l
（2H）]
4步Runge-Kutta法步长示意图
步1
n
步1x
R（
n,Cfn,Hn）
步2
n
步2x
R（
步1,Cf步1,H步1）
nrnfiff!
步m
n
步mx
R（
步m-1,Cf步m-1,H步m-1）
步
x
2^nax
m
解常微分方程的Runge-Kutta多步法:
（4）根据最后解得的附面层动量损失厚度0计算机翼上下表面的摩擦阻力。
DragDrag上Drag下
（5）利用整个计算分析系统，对不同设计方案的机翼开展摩擦阻力的对比分析。
由计算得到的形状因子说明各个方案气流分离情况（以H>3.55为标准）。
三、设计程序
functionOUTS=Drag_Airfoil
wave
%%%GenerictimemarchingcodesolvingthePDEforonedimensional
%%%WrittenbyHuangGuoping,2008/5/4
nmax=19;
%inputthedataofanairfoil
[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';
%plottheshapeofairfoil
plotfoil（XU,YU,XL,YL）;
%computetheboundarylayerofairfoil'suppersurfacelengthU
（1）=0;thetaU
（1）=0;CfU
（1）=0;HU
（1）=1;
forn=2:
N_U
dx（n）=dis（XU,YU,n）;
lengthU（n）=lengthU（n-1）+dx（n）;
ifn==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
%computetheboundarylayerofairfoil'slowersurface
lengthL
（1）=0;thetaL
（1）=0;CfL
（1）=O;HL
（1）=1;
forn=2:
N_L
dx（n）=dis（XU,YU,n）;
lengthL（n）=lengthL（n-1）+dx（n）;
ifn==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
%plottheresultsofairfoil
%outputtheUppersurface'sparameters
plotResults（lengthU,VU/Vsound,thetaU/（Chord*0.01）,HU）;
%outputtheLowersurface'sparameters
plotResults（lengthL,VL/Vsound,thetaL/（Chord*0.01）,HL）;
4/17
%ObtainthefrictionalDarg
DragU=thetaU（N_U）*Span*Density*Vupstream*Vupstream;
DragL=thetaL（N_L）*Span*Density*Vupstream*Vupstream;
Drag=DragU+DragL;
%ENDOFMAIN
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function
[Density,Tem,Vupstream,Chord,Span,DataU,N_U,DataL,N_L]=inputData（）
%N=input（'enternoofgridpoints__'）;
file1=fopen（
'foil-0.dat','r'）;
ccc=fscanf（file1.
'%7f%7f%7f%7f%7f',[51]）';
Density=ccc
（1）;Tem=ccc
（2）;Vupstream=ccc（3）;Chord=ccc（4）;
Span=ccc（5）;
tempc=fscanf（file1,
'%25c',[11]）;
N_U=fscanf（file1.
'%5i',[11]）;
tempc=fscanf（file1,
'%25c',[11]）;
N_L=fscanf（file1,
'%5i',[11]）;
fclose（file1）;
file2=fopen（'foil-U.dat','r'）;
DataU=fscanf（file2,'%8f%8f%7f%6f',[4N_U]）';
fclose（file2）;
file3=fopen（'foil-L.dat','r'）;
[4N_L]）';
DataL=fscanf（file3,'%8f%8f%7f%6f
fclose（file3）;
%END
functionmiu=Sutherland（Tem）
miu0=1.4587e-6;Tem0=110.4;
miu=miu0*（（Tem）A1.5）/（Tem+Tem0）;
%END
functiondistance=dis（X,Y,n）
distance=sqrt（（X（n）-X（n-1））A2+（Y（n）-Y（n-1）F2）;
%END
functionplotfoil（XU,YU,XL,YL）
figure
holdon;
plot（XU,YU,'-o'）;
plot（XL,YL,'-o'）；
axis'equal'
holdoff
%END
functionplotResults（L,V,theta,H）
figure
plot（L,V,'-d'）;
figure
holdon;
plot（L,theta,'-d'）;
plot（L,H,'-o'）;
holdoff;
%END
function[theta,H]=BoundaryLayer_Flatplate（length,V,Density,miu）
Rel=Density*V*length/miu;
%BlasuisSolutionforlaminarflow
theta=0.664*length/sqrt（Rel）;
Cf=0.664/sqrt（Rel）;
H=2.59;
%{
%Algorithmforturbulentflow
theta=0.0142*（RelA（6/7））*miu/（Density*V）;
Cf=0.026*（RelA（-1/7））;
H=1.375;
%}
%END
function[LL,HH]=Thwaites（Lamda）
ifLamda>0.25
Lamda=0.25;
else
ifLamda<=-0.09;
Lamda=-0.09;
end
end
ifLamda>=0
LL=0.22+1.57*Lamda-1.8*LamdaA2;
HH=2.61-3.75*Lamda-5.24*LamdaA2;
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）
%SolutionofRunge-Kuttamethod
Sita
（1）=theta1;
CF=O;
HH=1;
Sita
（2）=Sita
（1）+dx*（CF/2-（2+HH）*Sita
（1）*（V（n）-V（n-1））/V（n）/dx）/8;fori=2:
1:
4
Lamda=Density*（Sita（i）A2）*（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
A（4-i））；
end
H=HH;
Cf=CF;
theta=Sita（i）;
%END
function[theta2,HH]=
BoundaryLayerEquation1（dx,n,V,Density,miu,theta1）
%SolutionofEulermethodcoupledwithIteration
theta2=thetal;theta_error=theta1;
Iter=0;
while（theta_error>abs（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丄amda））;
[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
四、数据分析
1、弯度不变，弦长和翼展变化：
（1）第一组数据：
Camber=3.0%chord；
Thickness=9.0%chord；
Span=0.5m。
Chord=0.3m；
图中figure1-5分别为：
a）机翼外形；
b）上表面的速度分布；
c）下表面的速度分布；
d）上表面H因子和Theta的变化;
e）下表面H因子和Theta的变化。
figure4
Siiifi■■兰ZJlit.Ui£>tUm♦霸凹
详细输出数据:
lJd・k
kDW-dBQ■口
上表面
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
out=7.0000
7.6396
0.0001
0.0007
2.8416
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
附面层出现分离
F表面:
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
out=5.0000
4.6670
0.0000
0.0035
2.4251
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
附面层出现分离
DFiguie1
轴旧«wei和旳爲a®jajD直回eib网m的
JUi命=、£⑥©帚-Q□Q
DragU=3.8960
DragL=1.8175
Drag=5.7134
分析：
以上数据表明，所设计机翼其上表面在第11个点就出现附面层分离，下表面只有最
后两点出现附面层分离，因此所设计机翼需要进行改进。
改进情况见后三组数据。
（2）第二组数据：
Camber=3.0%chord；
Thickness=9.0%chord；
Chord=0.355m；
Span=0.455m。
图中figure1-5分别为：
a）机翼外形；
b）上表面的速度分布；
c）下表面的速度分布；
d）上表面H因子和Theta的变化;
e）下表面H因子和Theta的变化。
okrir|jj口up血
QK«1«1<5a.2。
站mF93&
Drag=5.3640
（3）结论：
比较上图可知，当camber不变的时候，chord和span的变化对分离的点的影响较小，基本可以忽略不计。
空气阻力却与chord和span有着密切的关系，当chord增大,span减小时，阻力会随之减小，反之则阻力随之增大。
具体原因是由于机翼前缘附面层较薄，因此速度梯度较大，所以机翼前缘的粘性阻力较大，机翼沿流线方向向后则空气阻力随之减小。
因此，机翼弦长较短，翼展较大时，相对的机翼前缘就比较长，所以空气阻力就较大，反之则空气阻力较小。
2、弯度变化，弦长和翼展变化忽略时：
（1）第一组数据：
Camber=1.0%chord；
Thickness=8.0%chord；
Chord=0.355m；
Span=0.7m。
图中figure1-5分别为：
a）机翼外形；
b）上表面的速度分布；
c）下表面的速度分布；
d）上表面H因子和Theta的变化;
e）下表面H因子和Theta的变化。
Drag=1.6908
（2）第二组数据：
Camber=5.0%chord；
Thickness=8.0%chord；
Chord=0.250m;
Span=0.500m。
Drag=2.8527
（3）结论：
由上述过程中，可以发现当camber增大，Thickness相对不变时，其上表面的分离位置向后移动，但是下表面的中部会出现分离点，因此结论是设计机翼时需要根据具体情况设计，并且需要进行多次的反复的修改和优化，以达到最优的设计。
camber
增大，其摩擦阻力成增大趋势；根据第一二组数据可知：
当chord增大，span减小
时，其摩擦阻力成减小趋势。
五、心得体会
经过将近一星期的努力我终于独立完成了这次粘性流体力学的大作业！
相信一开
始大家都一样毫无头绪，但认真地看了一遍老师的大作业讲解PPT后，才发现其实所
有步骤老师都用图文精心地准备好了。
知道了要做什么，我首先又恶补了一番matlab，
把老师给的源程序一运行发现错误甚多。
认真检查了才发现少了好多end,函数都没
结束•。
之后在老师留“？
”的地方补好函数后运行就得到了5张图。
在分析的过程中，
我觉得自己对粘性流体力学这门课学地更加透彻了，同时也对老师留的Foil.html工具
产生了浓厚的兴趣，甚至登陆了NASA官网逛了一圈，收获颇丰。
相信这次的大作业
给大家都留下了深刻的印象！
升级会员 