计算流体力学课程大作业.docx
- 文档编号:7199171
- 上传时间:2023-01-21
- 格式:DOCX
- 页数:15
- 大小:802.12KB
计算流体力学课程大作业.docx
《计算流体力学课程大作业.docx》由会员分享,可在线阅读,更多相关《计算流体力学课程大作业.docx(15页珍藏版)》请在冰豆网上搜索。
计算流体力学课程大作业
《计算流体力学》课程大作业
基于涡量-流函数法的不可压缩方腔驱动流问题数值模拟
张伊哲航博101
1、引言和综述
2、问题的提出,怎样使用涡量-流函数方法建立差分格式
3、程序说明
4、计算结果和讨论
5、结论
1引言
虽然不可压缩流动的控制方程从形式上看更为简单,但实际上,目前不可压缩流动的数
值方法远远不如可压缩流动的数值方法成熟。
考虑不可压缩流动的N-S方程:
r—0
—'(UU)=f-丄\PU(1.1)
其中、是运动粘性系数,认为是常数。
将方程组写成无量纲的形式:
[U=0
:
U1(1.2)
—(UU)=f」P丄U
ftRe
其中Re是雷诺数。
从数学角度看,不可压缩流动的控制方程中不含有密度对时间的偏导数项,方程表现出
椭圆-抛物组合型的特点;从物理意义上看,在不可压缩流动中,压力这一物理量的波动具有无穷大的传播速度,它瞬间传遍全场,以使不可压缩条件在任何时间、任何位置满足,这
就是椭圆型方程的物理意义。
这就造成不可压缩的N-S方程不能使用比较成熟的发展.型.偏
微分方程的数值求解理论和方法。
如果将动量方程和连续性方程完全耦合求解,即使使用显示的离散格式,也将会得到一
个刚性很强的、庞大的稀疏线性方程组,计算量巨大,更重要的问题是不易收敛。
因此,实际应用中,通常都必须将连续方程和动量方程在一定程度上解耦。
目前,求解不可压缩流动的方法主要有涡量-流函数法,SIMPLE法及其衍生的改进方
法,有限元法,谱方法等,这些方法各有优缺点。
其中涡量-流函数法是解决二维不可压缩
流动的有效方法。
作者本学期学习了研究生计算流体课程,为了熟悉计算流体的基本方法,选择使用涡量-流函数法计算不可压缩方腔驱动流问题,并且对于不同雷诺数下的解进行比较和分析,得出一些结论。
本文接下来的内容安排为:
第2节提出不可压缩方腔驱动流问题,并分析该问题怎样使
用涡量-流函数方法建立差分格式、选择边界条件。
第3节介绍程序的结构。
第4节对于不
同雷诺数下的计算结果进行分析,并且与UGHIA等人【1】的经典结论进行对比,评述本
文所采用的计算方法。
第五节给出结论。
2问题的提出和分析
2.1经典方腔驱动流问题
考虑如下图所示的长度为1的正方形腔体,腔体上有一平板以速度U=1运动,其它三
边为固壁条件。
图1•方腔驱动流示意图
顶盖方腔驱动流问题是个很经典的问题,常常用于验证不可压缩流动数值方法的正确
性。
U.GHIA等人于1982年发表的一篇文献(见文献【1】)计算了Re从100到104的流动
结果,其结果得到广泛的认同。
2.2涡量-流函数方法简介
涡量-流函数法的基本思想是引入涡量「和流函数•:
引入涡量,可以消去方程中的压力项,而引入流函数,可以使连续方程自然满足。
下面对该方法进行简单推导:
(1.2)写成分量形式:
考虑二维问题,将式
(1-3)
-U:
vc
0
.:
t亍
可见,连续性方程(1.3)自然成立。
’-:
与•,的关系为
(1.8)
式(1.6)~(1.8)构成了一个封闭的方程组,由(1.6)计算出涡量,再由(1.8)式计算出流函数,
利用(1.7)式计算出速度。
这个方程组的特点是求解速度的时候完全不用考虑压力项。
若还需
要求解压力场,则可以把式(1.4)对X求偏导数,式式(1.5)对y求偏导数,二者求和后整理得到关于压力的Poisson方程
:
y;x
0.8
以上推导出的涡量
-流函数法在计算二维问题时很成功,但是三维流动的流函数没有直
观的物理意义,无法像二维流动一样直接定义,需要引入多个流函数,相应解多个Poisson
方程,计算量很大,并不实用。
对于本文的二维问题,该方法就简单易行。
2.3建立差分格式
2.3.1划分网格
方腔驱动流的流动区域很简单,均匀划分为正方形的结构网格即可,存储网格时,x方
向使用标号i表示,y方向使用标号j表示,x和y方向的最大网格点标号分别为M和N。
对于Re小于等于1000的情况,使用100*100网格,Re大于1000后的情况,使用256*256网格。
计算域如图2所示:
■□■!
5a!
!
5O!
!
SH3:
!
S!
!
OS!
!
2E2E:
:
!
!
SO=!
SO!
!
S!
'E!
S!
OOO3:
2O!
!
5O!
!
■:
E:
・bi・・・■・・・■・・・・bi・・・■■・・・■■・・■■・・・im・・・r・・■・・■・・・m・・・■■・・・・■■・・・■■・・・■■・・・・an・・・■・・・■■・・■■・・・ta
anSi3aE!
aQEI13EiSiHiEiaDEi33Eis:
HiEiHisEDIliDEiEUEiESDiiiDEiiDEia3EiSi3DBiaDi:
UEiSi3
■□■!
S!
3!
!
SO!
!
2:
:
!
S!
!
OS!
H3!
!
5!
!
E!
!
SOS!
SO!
!
S!
2!
!
S!
O!
!
:
nE5!
[5S!
!
SE!
!
22E
iEiiiiMiiiMiiiiiiiiiiiiiiEiiiniiiniiiiiiiiiMiiiaiiiniiiiMiiiMiiifiiiiiiiiiiMiiiMiiuiiiii
aM;:
3u;:
3O;:
33E:
:
:
3ZiS:
3OE!
33;ii:
jZi:
:
4L;iaOSi3u;:
S:
j»:
ain;;3ui;:
iO;Iia;:
:
:
au;:
SnE:
3jEii:
j
图2.100*100的均分网格
232建立差分方程
由于本题关注的是方腔内部的流动状态,对于压力分布不关心,因此不用建立压力的差
分方程。
涡量的对流扩散方程(1.6)使用FTCS格式离散得到:
该差分格式时间方向为1阶精度,空间方向为2阶精度。
在(1.10)中,速度分量取的是n时刻的值,已经对方程进行了线性化处理。
流函数的Poisson方程中,二阶导数都用中心差分离散:
这种中心差分可达到二阶精度。
2.3.3设定边界条件
(1)速度和流函数的边界条件
由于沿着壁面是一条流线,所以流函数在边界是常值,可以取为o;速度在边界满足无
滑移条件。
上边界(i=O~M,j=N):
Ui,N=U=1,v,n=0;Vj,N=0
下边界(i=0〜M,j=0):
5,0=0,w,0=0;屮i,0=0
左边界(j=1〜N—1,i=0):
u,0j=0qVj0;即0,j=0
右边界(j=1〜N-1,i=M):
Um,j=0,Vm,j=0;M,j=0
(2)涡量的边界条件
'叩「2:
:
冋
左边界(j=1〜N—1,i=0):
co0j二一,这里引入了虚拟网格点(-1,j),
4
注意到1;;j,所以訥0」=红
j七〜'x
2屮
同理可得,右边界(j=1〜N—1,i=M):
⑷二一字
Z
下边界(i=0〜M,j=0):
偽0二生
2U,n」2逍
下面考察构造的这种涡量边界条件的精度:
比如对于下边界,流函数的Taylor展开为
对于其他边界的精度推导是类似的。
式简单,还能有2阶精度。
3编程计算
3.1程序的结构
程序流程图如下图所示:
图3•流程图
3.2程序说明
示法,一步就可以将•打;1求出。
(1.12)
2、对于流函数场的计算:
流函数满足泊松方程,差分形式为:
n1_2n1.,_:
n1.n1_2[-■nd..n1
i1,j_2i,ji4,ji,jd_2i,ji,j4n1
"T=O..
Ax2Ay2i,j
求解时要使用JACOBI迭代,由于歆》:
:
y,可以写成
■7.n1,(k1)=丄tn1,(k)胪:
yn1,(k)胪:
yn1,(k):
:
冋n1,(k)_山2.n1
i,j-.i1,jiJ,ji,j1i,j-1-h'i,j
4
其中(k+1)表示第k+1次迭代,当|旷+1,(k+1)-旷皿丄<E时,认为精度达到,退出迭代。
将此时的»n+1,(k+1)作为n+1时刻的流函数场»n+1。
3、速度场的求解
6屮云屮
由于v,U,所以
:
x:
y
4结果分析
图4是Re=100,400,1000,3200,5000,7500,10000时的流线图,每个雷诺数中,上图是本文计算的图,下图是文献1中的结果。
与文献中图形进行比较,流线图很相似,计算结果是可信的。
从图4中看出,方腔驱动流中,Re=100,400,1000时,只在左右两个下角落出现二次
涡,Re=3200时,在靠近顶盖(注意不是在顶盖上)的左边壁面又出现了一个二次涡,Re=7500
时,右边的下角落存在两个二次涡,所有二次涡的强度都随着Re增大而增大。
中央主涡的
中心当Re=100时偏向右上方,随着Re增大逐渐移动到方腔的几何中心上,当Re=3200也
就是当左边壁面出现二次涡之后,主涡的中心几乎保持不变。
可以预测,当Re接着增大,
在现有二次涡的地方会出现更多的二次涡。
图4.不同Re时的流线图
图5显示了不同Re数中,在通过方腔几何中心的竖直线上速度u的分布情况,从图5
中看出,随着Re增大,边界层越来越薄,Re>5000后,边界层变薄的速率就很小了。
高Re
数时,方腔中部的速度分布几乎是线性的。
Re>3200时,在靠近y=1处,u的分布会出现一
个小凸起,这个现象也被之前的文献提到过。
图5.不同Re时过方腔几何中心的竖直线上u的分布
图6显示了Re=100,400,1000,5000,10000时,文献1中计算出的竖直中心线上u的
分布和本文计算的结果,其中,Re=100,400,1000时,本文使用100*100均分网格,文献1使用129*129优化网格;Re=5000和10000时,本文使用257*257均分网格,文献1使用257*257优化网格。
可见,当Re比较低的时候,本文的计算是很准确的;当Re增加至5000,
即使将网格增加至257*257,仍然存在较大误差,尤其是在靠近y=1顶盖处,因为本文的网
格是均分的,即使加密,效果也远远不如文献1中的优化网格。
-U/l-U20U.dU.b1
u
高Re时,计算时间较长,如下表所示。
Re
00
1
400
1000
3200
5000
750
10
0
000
CPU
间(s)
时
0
9
166
208
810
940
5
110
09
15
5结论
本文通过涡量-流函数法对不同雷诺数下的方腔驱动流进行了模拟,得到结论如下:
对于方腔驱动流:
1、随着Re增加,二次涡逐个出现。
2、随着Re增加,主涡的中心逐渐移动到方腔的几何中心,Re>3200后,主涡的位置几乎不变。
3、随着Re增加,边界层越来越薄。
高Re数时,方腔中部的速度分布几乎是线性的,在靠近y=1处,u的分布会出现一个小凸起。
对于本文使用的算法
1、Re<3200时,计算结果精良,耗时短。
2、由于网格是均分,高Re时计算结果在边界处有较大误差。
3、可以改进计算流函数泊松方程时的方法,比如使用共轭梯度法,使得收敛加快,节省计算时间。
参考文献:
[1]GHIAU,GHIAKN,SHINCT.High-ResolutionsforincompressibleflowusingtheNavier-StokesEquationsandamultigridmethod]J].JComputPhys,1982,48:
3872411.
后记:
本学期学习了李嵩老师开设的计算流体课程,受益颇多。
对于计算流体的基本方法和思路有所了解,于是就编写了简单的程序作为大作业。
由于第一
次大作业的选题过于难,没有做出,跟李老师沟通后,老师非常理解,宽限我一些时间,在此衷心感激老师的理解和宽容,就老师一学期的谆谆教诲表示感谢。
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 计算 流体力学 课程 作业
![提示](https://static.bdocx.com/images/bang_tan.gif)