计算流体力学课程大作业
- 格式:doc
- 大小:1.08 MB
- 文档页数:12
《计算流体力学》课程大作业
——基于涡量-流函数法的不可压缩方腔驱动流问题数值模拟
张伊哲 航博101
1、 引言和综述
2、 问题的提出,怎样使用涡量-流函数方法建立差分格式
3、 程序说明
4、 计算结果和讨论
5、 结论
1引言
虽然不可压缩流动的控制方程从形式上看更为简单,但实际上,目前不可压缩流动的数值方法远远不如可压缩流动的数值方法成熟。
考虑不可压缩流动的N-S 方程:
01()P t νρ∇⋅=⎧
⎪
∂⎨+∇⋅=-∇+∆⎪∂⎩
U U
UU f U (1.1)
其中ν是运动粘性系数,认为是常数。将方程组写成无量纲的形式:
01()Re P t
∇⋅=⎧⎪
∂⎨+∇⋅=-∇+∆⎪∂⎩U U
UU f U (1.2) 其中Re 是雷诺数。 从数学角度看,不可压缩流动的控制方程中不含有密度对时间的偏导数项,方程表现出椭圆-抛物组合型的特点;从物理意义上看,在不可压缩流动中,压力这一物理量的波动具有无穷大的传播速度,它瞬间传遍全场,以使不可压缩条件在任何时间、任何位置满足,这就是椭圆型方程的物理意义。这就造成不可压缩的N-S 方程不能使用比较成熟的发展型...偏微分方程的数值求解理论和方法。
如果将动量方程和连续性方程完全耦合求解,即使使用显示的离散格式,也将会得到一个刚性很强的、庞大的稀疏线性方程组,计算量巨大,更重要的问题是不易收敛。因此,实际应用中,通常都必须将连续方程和动量方程在一定程度上解耦。
目前,求解不可压缩流动的方法主要有涡量-流函数法,SIMPLE 法及其衍生的改进方法,有限元法,谱方法等,这些方法各有优缺点。其中涡量-流函数法是解决二维不可压缩流动的有效方法。作者本学期学习了研究生计算流体课程,为了熟悉计算流体的基本方法,选择使用涡量-流函数法计算不可压缩方腔驱动流问题,并且对于不同雷诺数下的解进行比较和分析,得出一些结论。
本文接下来的内容安排为:第2节提出不可压缩方腔驱动流问题,并分析该问题怎样使用涡量-流函数方法建立差分格式、选择边界条件。第3节介绍程序的结构。第4节对于不同雷诺数下的计算结果进行分析,并且与U.GHIA 等人【1】的经典结论进行对比,评述本
文所采用的计算方法。第五节给出结论。
2问题的提出和分析
2.1经典方腔驱动流问题
考虑如下图所示的长度为1的正方形腔体,腔体上有一平板以速度U=1运动,其它三边为固壁条件。
图1.方腔驱动流示意图
顶盖方腔驱动流问题是个很经典的问题,常常用于验证不可压缩流动数值方法的正确性。U.GHIA 等人于1982年发表的一篇文献(见文献【1】)计算了Re 从100到410的流动结果,其结果得到广泛的认同。
2.2涡量-流函数方法简介
涡量-流函数法的基本思想是引入涡量ω和流函数ψ:引入涡量,可以消去方程中的压力项,而引入流函数,可以使连续方程自然满足。下面对该方法进行简单推导:
考虑二维问题,将式(1.2)写成分量形式:
222222220 (1.3)1 (1.4)Re 1 (1.5)Re u v
t t u u u P u u u v t x y x x y v v v P v v u v t x y y x y ∂∂+=∂∂⎛⎫∂∂∂∂∂∂++=-++ ⎪∂∂∂∂∂∂⎝⎭⎛⎫∂∂∂∂∂∂++=-++ ⎪∂∂∂∂∂∂⎝⎭
式(1.4)对y 求偏导数减去式(1.5)对x 求偏导数,考虑到=
u v
y x
ω∂∂-∂∂,推导出涡量满足
的方程为22221Re u v t x y x y ωωωωω⎛⎫∂∂∂∂∂++=+ ⎪∂∂∂∂∂⎝⎭
(1.6) 然后引入流函数ψ,定义为
,v u x y
ψψ
∂∂=-=∂∂ (1.7) 可见,连续性方程(1.3)自然成立。ψ与ω的关系为
222
2x y
ψψ
ω∂∂+=∂∂ (1.8) 式(1.6)~(1.8)构成了一个封闭的方程组,由(1.6)计算出涡量,再由(1.8)式计算出流函数,利用(1.7)式计算出速度。这个方程组的特点是求解速度的时候完全不用考虑压力项。若还需要求解压力场,则可以把式(1.4)对x 求偏导数,式式(1.5)对y 求偏导数,二者求和后整理得到关于压力的Poisson 方程
22
2
2u u v v P x y x y ⎧⎫⎛⎫∂∂∂∂⎪⎪⎛⎫∇=-++⎨⎬ ⎪ ⎪∂∂∂∂⎝⎭⎝⎭⎪⎪⎩
⎭ (1.9)
以上推导出的涡量-流函数法在计算二维问题时很成功,但是三维流动的流函数没有直
观的物理意义,无法像二维流动一样直接定义,需要引入多个流函数,相应解多个Poisson 方程,计算量很大,并不实用。对于本文的二维问题,该方法就简单易行。
2.3建立差分格式
2.3.1划分网格
方腔驱动流的流动区域很简单,均匀划分为正方形的结构网格即可,存储网格时,x 方向使用标号i 表示,y 方向使用标号j 表示,x 和y 方向的最大网格点标号分别为M 和N 。对于Re 小于等于1000的情况,使用100*100网格,Re 大于1000后的情况,使用256*256网格。计算域如图2所示:
图2.100*100的均分网格
2.3.2建立差分方程
由于本题关注的是方腔内部的流动状态,对于压力分布不关心,因此不用建立压力的差分方程。
涡量的对流扩散方程(1.6)使用FTCS 格式离散得到:
1,,1,1,,1,1
,,1,,1,,1
,,1
2
2
22221()
Re
()
()
n n n n n n i j i j
i j i j
i j i j n n i j
i j
n n n n n n i j
i j i j
i j i j i j u
v
t x
y
x y ωωωωωωωωωω
ωω++-+-+-+----++∆∆∆-+-+=∆∆ (1.10)
该差分格式时间方向为1阶精度,空间方向为2阶精度。在(1.10)中,速度分量取的是n 时刻的值,已经对方程进行了线性化处理。
流函数的Poisson 方程中,二阶导数都用中心差分离散:
1111111,,1,,1,,1
1
,2
2
22n n n n n n i j i j i j i j i j i j n i j x y ψψψψψψω+++++++-+-+-+-++
=∆∆ (1.11)
这种中心差分可达到二阶精度。
2.3.3设定边界条件
(1)速度和流函数的边界条件
由于沿着壁面是一条流线,所以流函数在边界是常值,可以取为0;速度在边界满足无滑移条件。
上边界(0~,i M j N ==): ,,1,0i N i N u U v ===;,0i N ψ= 下边界(0~,0i M j ==): ,0,00,0i i u v ==;,00i ψ= 左边界(1~1,0j N i =-=):0,0,0,0j j u v ==;0,0j ψ= 右边界(1~1,j N i M =-=):,,0,0M j M j u v ==;,0M j ψ= (2)涡量的边界条件
根据涡量的定义=
u v y x ω∂∂-∂∂,在上下边界,0v
x ∂=∂,所以22=u y y ψω∂∂=∂∂;在左右边界,0u
y
∂=∂,所以22=v x x ψω∂∂-=∂∂。;
左边界(1~1,0j N i =-=):
1,0,1,0,2
2=,j j j
j x ψψψω--+∆这里引入了虚拟网格点(-1,j ),