南京航空航天大学.docx
- 文档编号:491024
- 上传时间:2022-10-10
- 格式:DOCX
- 页数:29
- 大小:185.25KB
南京航空航天大学.docx
《南京航空航天大学.docx》由会员分享,可在线阅读,更多相关《南京航空航天大学.docx(29页珍藏版)》请在冰豆网上搜索。
南京航空航天大学
南京航空航天大学
计算空气动力学程序报告
姓名:
戴秋敏
学号:
S0701362
专业:
人机与环境工程
航空宇航学院
2008-01-17
由下图所示的翼型模型,根据已给定的网格,并编写程序计算在一定条件下机翼的升力系数和阻力系数。
图1:
翼型模型
本程序使用fortran语言,机翼网格采用最简单的C结构网格,其特点是程序编写简单,网格节点有很好的连贯性。
离散方法则采用格心格式,而对流通量的计算采用教材中素哟介绍的最基本的AUSM格式,未采用AUSM改进格式,界面左右变量的插值采用守恒变量。
网格边界处采用两层虚拟网格,虚拟单元中:
远场取来流值,物面采用对称处理,cutline处直接把相对应的单元值传递过去。
因为未在计算过程中用最小时间步长而不是当地时间步长,并且未采取各加速收敛措施,所以程序的收敛性不是很好,达不到很高的精度。
如下所示的已给定的机翼模型网格:
图2:
翼型网格
计算流体力学基本原理
守恒形式的Euler方程:
(1)
其中x和y代表笛卡儿坐标系。
W是守恒变量。
(2)
F,G表示通量
(3)
P,H和E表示密度,压强,单元总焓和单元总能量。
U,V表示笛卡儿坐标系下的速度矢量。
这些量由理想气体的单位体积的总能量和总焓相互联系。
(4)
(5)
程序计算结果分析
1.马赫数:
0.85,夹角:
2°,压力100000pa,温度:
288k.
密度和马赫数分布:
图6:
密度分布图7:
马赫数分布
压力和温度分布
图8:
压力分布图9:
温度分布
升力和阻力系数
升力系数:
1.654012693519855E-004
阻力系数:
7.913074813368969E-002
附:
C用AUSM格式求解二维Euler方程
IMPLICITDOUBLEPRECISION(A-H,O-Z)!
I-N规则
PARAMETER(IX=210,IY=50)!
宏定义值要比193,33大
DOUBLEPRECISIONMA
DIMENSIONRK(4)!
荣格-库塔
COMMON/INPUTD/MA,AFA,PIN,TIN,CFL,EPSL,ICYMAX,IREAD
C马赫数,迎角,来流压力,来流温度,cfl值,收敛标准,最大迭代数,是否接着以前的计算
COMMON/AIR/GAMA,R
C空气比热比,气体常数
COMMON/MESH/NDI,NDJ,NCI,NCJ,IS1,IS2
ci,j方向的网格节点数;i,j方向的单元序号最大值;i,j方向物面开始,结束网格点序号
COMMON/GRID/X(IX,IY),Y(IX,IY),X0(IX,IY),Y0(IX,IY)
c网格点坐标,单元中心点坐标
COMMON/VECTOR/SIU(IX,IY,2),SJU(IX,IY,2),SI(IX,IY),SJ(IX,IY),
&SCI(IX,IY,2),SCJ(IX,IY,2),VOL(IX,IY)
ci,j方向界面单位法矢、面积,单元中心i,j方向边单位法矢,单元面积
COMMON/VARI/CV(IX,IY,4),DV(IX,IY,3),W(IX,IY,4),CVIN(4),
&VIN,ROUIN,PIN0,ROUIN0,UINF0,VIN0,AIN0
c守恒变量,原始变量(p,t,a),"inf"表示来流,"0"表示无量刚值
COMMON/TS/TSTEP(IX,IY),TTSTEP
COMMON/RE/RES(IX,IY,4),EPSM(10000)
DATARK/0.0833D0,0.2069D0,0.4265D0,1.D0/!
每步龙格库塔系数
WRITE(*,*)'SOLVE2-DEULEREQUATION'
WRITE(*,*)'PRESSENTERTOSTART'
READ(*,*)
CALLINPUT!
读输入文件
CALLINITIAL!
流场初始化
CALLGRIDD!
计算网格几何参数
CALLBOUNDARY!
设置边界条件
IF(IREAD.EQ.1)THEN!
判断是接着算还是从头开始算
OPEN(5,FILE='temp.dat')
READ(5,*)(((CV(I,J,K),K=1,4),I=1,NCI),J=1,NCJ)
CLOSE(5)
ENDIF
ISTEP=0
10ISTEP=ISTEP+1
CALLTIMESTEP!
时间步长
DOIRK=1,4!
四步R-K推进
CALLRESID!
求残值
DOJ=3,NCJ
DOI=3,NCI
DOK=1,4
CV(I,J,K)=W(I,J,K)-RK(IRK)*TTSTEP*RES(I,J,K)/VOL(I,J)!
流场解更新
ENDDO
ENDDO
ENDDO
DOJ=3,NCJ
DOI=3,NCI
W1=CV(I,J,1)
W2=CV(I,J,2)
W3=CV(I,J,3)
W4=CV(I,J,4)
CALLTRANS2(W1,W2,W3,W4,DV(I,J,1),DV(I,J,2),DV(I,J,3))
ENDDO
ENDDO
CALLBOUNDARY!
设置边界
ENDDO
EPSM(ISTEP)=0.
DOJ=3,NCJ
DOI=3,NCI
EPSM(ISTEP)=EPSM(ISTEP)+(CV(I,J,1)-W(I,J,1))**2!
求解的变化大小用来判断收敛
ENDDO
ENDDO
W=CV
IF(ISTEP==1)EPS1=EPSM(ISTEP)
EPSM(ISTEP)=SQRT((EPSM(ISTEP)+1.E-32)/EPS1)
WRITE(*,*)ISTEP,EPSM(ISTEP)
IF(EPSM(ISTEP).GT.EPSL.AND.istep.LT.ICYMAX)THEN
GOTO10
ELSE
CALLOUTPUT(ISTEP)
ENDIF
END
C----------------------------------------------------------------
SUBROUTINEINPUT!
输入
IMPLICITDOUBLEPRECISION(A-H,O-Z)
PARAMETER(IX=210,IY=50)
DOUBLEPRECISIONMA
DIMENSIONRK(4)
COMMON/INPUTD/MA,AFA,PIN,TIN,CFL,EPSL,ICYMAX,IREAD
C马赫数,迎角,来流压力,来流温度,cfl值,收敛标准,最大迭代数,是否接着以前的计算
COMMON/AIR/GAMA,R
C空气比热比,气体常数
COMMON/MESH/NDI,NDJ,NCI,NCJ,IS1,IS2
ci,j方向的网格节点数;i,j方向的单元序号最大值;i,j方向物面开始,结束网格点序号
COMMON/GRID/X(IX,IY),Y(IX,IY),X0(IX,IY),Y0(IX,IY)
c网格点坐标,单元中心点坐标
COMMON/VECTOR/SIU(IX,IY,2),SJU(IX,IY,2),SI(IX,IY),SJ(IX,IY),
&SCI(IX,IY,2),SCJ(IX,IY,2),VOL(IX,IY)
ci,j方向界面单位法矢、面积,单元中心i,j方向边单位法矢,单元面积
COMMON/VARI/CV(IX,IY,4),DV(IX,IY,3),W(IX,IY,4),CVIN(4),
&VIN,ROUIN,PIN0,ROUIN0,UINF0,VIN0,AIN0
c守恒变量,原始变量(p,t,a),"in"表示来流,"0"表示无量刚值
COMMON/TS/TSTEP(IX,IY),TTSTEP
COMMON/RE/RES(IX,IY,4),EPSM(10000)
OPEN(1,FILE='INPUT.DAT')
c马赫数,迎角,比热比,来流压力,温度,CFL数,收敛标准,最大迭代数,是否接着以前的解算
READ(1,*)MA!
马赫数
READ(1,*)AFA!
迎角
READ(1,*)GAMA!
比热比
READ(1,*)PIN!
来流压力
READ(1,*)TIN!
温度
READ(1,*)CFL!
cfl数
READ(1,*)EPSL!
收敛标准
READ(1,*)ICYMAX!
最大迭代数
READ(1,*)IREAD!
是否接着以前的解算
CLOSE
(1)
OPEN(2,FILE='NACA0012.DAT')
READ(2,*)NDI,NDJ!
i,j方向网格点数193,33
READ(2,*)IS1,IS2!
物面开始,结束网格点序号33,161
READ(2,*)((X(I,J),Y(I,J),I=1,NDI),J=1,NDJ)!
先执行I,再执行J
C(1,1),(2,1),(3,1)......
CLOSE
(2)
END
C------------------------------------------------------------------
SUBROUTINEINITIAL!
流场初始化
IMPLICITDOUBLEPRECISION(A-H,O-Z)
PARAMETER(IX=210,IY=50)
DOUBLEPRECISIONMA
DIMENSIONRK(4)
COMMON/INPUTD/MA,AFA,PIN,TIN,CFL,EPSL,ICYMAX,IREAD
C马赫数,迎角,来流压力,来流温度,cfl值,收敛标准,最大迭代数,是否接着以前的计算
COMMON/AIR/GAMA,R
C空气比热比,气体常数
COMMON/MESH/NDI,NDJ,NCI,NCJ,IS1,IS2
ci,j方向的网格节点数;i,j方向的单元序号最大值;i,j方向物面开始,结束网格点序号
COMMON/GRID/X(IX,IY),Y(IX,IY),X0(IX,IY),Y0(IX,IY)
c网格点坐标,单元中心点坐标
COMMON/VECTOR/SIU(IX,IY,2),SJU(IX,IY,2),SI(IX,IY),SJ(IX,IY),
&SCI(IX,IY,2),SCJ(IX,IY,2),VOL(IX,IY)
ci,j方向界面单位法矢、面积,单元中心i,j方向边单位法矢,单元面积
COMMON/VARI/CV(IX,IY,4),DV(IX,IY,3),W(IX,IY,4),CVIN(4),
&VIN,ROUIN,PIN0,ROUIN0,UINF0,VIN0,AIN0
c守恒变量,原始变量(p,t,a),"inf"表示来流,"0"表示无量刚值
COMMON/TS/TSTEP(IX,IY)
COMMON/RE/RES(IX,IY,4),EPSM(10000)
R=287.06D0!
气体
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 南京航空航天 大学