计算流体力学课程大作业

《计算流体力学》课程大作业——基于涡量-流函数法的不可压缩方腔驱动流问题数值模拟张伊哲 航博1011、 引言和综述2、 问题的提出，怎样使用涡量-流函数方法建立差分格式3、 程序说明4、 计算结果和讨论5、 结论1引言虽然不可压缩流动的控制方程从形式上看更为简单，但实际上，目前不可压缩流动的数值方法远远不如可压缩流动的数值方法成熟。

考虑不可压缩流动的N-S 方程：01()P t νρ∇⋅=⎧⎪∂⎨+∇⋅=-∇+∆⎪∂⎩U UUU f U (1.1)其中ν是运动粘性系数，认为是常数。

将方程组写成无量纲的形式：01()Re P t∇⋅=⎧⎪∂⎨+∇⋅=-∇+∆⎪∂⎩U UUU 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 vt 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 vy xω∂∂-∂∂，推导出涡量满足的方程为22221Re u v t x y x y ωωωωω⎛⎫∂∂∂∂∂++=+ ⎪∂∂∂∂∂⎝⎭(1.6) 然后引入流函数ψ，定义为,v u x yψψ∂∂=-=∂∂ (1.7) 可见，连续性方程(1.3)自然成立。

ψ与ω的关系为2222x yψψω∂∂+=∂∂ (1.8) 式(1.6)~(1.8)构成了一个封闭的方程组，由(1.6)计算出涡量，再由(1.8)式计算出流函数，利用(1.7)式计算出速度。

这个方程组的特点是求解速度的时候完全不用考虑压力项。

若还需要求解压力场，则可以把式(1.4)对x 求偏导数，式式(1.5)对y 求偏导数，二者求和后整理得到关于压力的Poisson 方程2222u 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,,12222221()Re()()n n n n n n i j i ji j i ji j i j n n i ji jn n n n n n i ji j i ji j i j i j uvt xyx y ωωωωωωωωωωωω++-+-+-+----++∆∆∆-+-+=∆∆ (1.10)该差分格式时间方向为1阶精度，空间方向为2阶精度。

在(1.10)中，速度分量取的是n 时刻的值，已经对方程进行了线性化处理。

流函数的Poisson 方程中，二阶导数都用中心差分离散：1111111,,1,,1,,11,2222n 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 ω∂∂-∂∂，在上下边界，0vx ∂=∂，所以22=u y y ψω∂∂=∂∂；在左右边界，0uy∂=∂，所以22=v x x ψω∂∂-=∂∂。

；左边界（1~1,0j N i =-=）：1,0,1,0,22=,j j jj x ψψψω--+∆这里引入了虚拟网格点（-1,j ），注意到0,1,1,0,002j j j j u x ψψψ-=⎧⎪-⎨==⎪∆⎩，所以1,0,22=j j x ψω∆。

同理可得，右边界（1~1,j N i M =-=）：1,,22=M j M j xψω-∆下边界（0~,0i M j ==）：,1,022=i i yψω∆上边界（0~,i M j N ==）：注意到,,1,1,012i N i N i N i Nu y ψψψ+-=⎧⎪-⎨==⎪∆⎩，所以,1,222=i N i N yyψω-+∆∆下面考察构造的这种涡量边界条件的精度：比如对于下边界，流函数的Taylor 展开为22334,1,023(,0)(,0)(,0)22334,1,023(,0)(,0)(,0)()2!3!()2!3!i i i i i i i i i i y y yO y yy y y y yO y y y y ψψψψψψψψψψ-∂∆∂∆∂=+∆+++∆∂∂∂∂∆∂∆∂=-∆+-+∆∂∂∂则42,1,0,1,1,0,12222(,0)2()2()i i i i i i i O y O y yy y ψψψψψψψ---++∆-+∂==+∆∂∆∆对于其他边界的精度推导是类似的。

所以构造的这种涡量边界条件很有优势——不仅形式简单，还能有2阶精度。

3编程计算3.1程序的结构程序流程图如下图所示：图3.流程图3.2程序说明时间步长选择0.001，足以满足稳定性条件。

1、对于涡量场的计算：根据差分方程1,,1,1,,1,11,,1,,1,,1,,22221()22Re ()()n n n n n n n n n n n ni j i ji j i ji j i j i j i j i j i j i j i j nn i ji juvtxy x y ωωωωωωωωωωωω++-+-+-+-----+-+++=∆∆∆∆∆该式中，速度分量取的是n 时刻的值，已经对方程进行了线性化处理。

所以直接使用显示法，一步就可以将1,n i j ω+求出。

1,,1,,1,,1,,11,1,,1,1,,22221()Re ()()22n n i j i j n n n n n n n n n n i j i j i j i j i j i j i j i j i j i j n n i j i j t u v x y x y ωωωωωωωωωωωω++-+-+-+-=+⎡⎤⎛⎫-+-+--∆-+⎢⎥ ⎪ ⎪∆∆∆∆⎢⎥⎝⎭⎣⎦（1.12）2、对于流函数场的计算：流函数满足泊松方程，差分形式为：1111111,,1,,1,,11,2222n n n n n n i j i j i j i j i j i j n i jx y ψψψψψψω+++++++-+-+-+-++=∆∆求解时要使用JACOBI 迭代，由于x y ∆=∆，可以写成()1,(1)1,()1,()1,()1,()21,1,1,,1,1,14n k n k n k n k n k n i ji j i j i j i j i j h ψψψψψω++++++++-+-=+++- 其中（k+1）表示第k+1次迭代，当1ε-<n+1,(k+1)n+1,(k)ψψ时，认为精度达到，退出迭代。

将此时的n+1,(k +1)ψ作为n+1时刻的流函数场n+1ψ。

3、速度场的求解 由于,v u x yψψ∂∂=-=∂∂，所以 ,1,11,1,,22i j i j i j i ju v yxψψψψ+-+---==-∆∆4结果分析图4是Re=100,400,1000,3200,5000,7500,10000时的流线图,每个雷诺数中，上图是本文计算的图，下图是文献1中的结果。