计算流体力学课程大作业2.pdf
计算流体力学课程大作业计算流体力学课程大作业基于涡量-流函数法的不可压缩方腔驱动流问题数值模拟张伊哲航博 1011、引言和综述2、问题的提出，怎样使用涡量-流函数方法建立差分格式3、程序说明4、计算结果和讨论5、结论1 1 引言引言虽然不可压缩流动的控制方程从形式上看更为简单，但实际上，目前不可压缩流动的数值方法远远不如可压缩流动的数值方法成熟。考虑不可压缩流动的 N-S 方程：U U 01U U(UUUU)f f PU Ut(1.1)其中是运动粘性系数，认为是常数。将方程组写成无量纲的形式：U U 0(1.2)1U U(UUUU)f f PU URe t其中 Re 是雷诺数。从数学角度看，不可压缩流动的控制方程中不含有密度对时间的偏导数项，方程表现出椭圆-抛物组合型的特点；从物理意义上看，在不可压缩流动中，压力这一物理量的波动具有无穷大的传播速度，它瞬间传遍全场，以使不可压缩条件在任何时间、任何位置满足，这就是椭圆型方程的物理意义。这就造成不可压缩的 N-S 方程不能使用比较成熟的发展型偏微分方程的数值求解理论和方法。如果将动量方程和连续性方程完全耦合求解，即使使用显示的离散格式，也将会得到一个刚性很强的、庞大的稀疏线性方程组，计算量巨大，更重要的问题是不易收敛。因此，实际应用中，通常都必须将连续方程和动量方程在一定程度上解耦。目前，求解不可压缩流动的方法主要有涡量-流函数法，SIMPLE 法及其衍生的改进方法，有限元法，谱方法等，这些方法各有优缺点。其中涡量-流函数法是解决二维不可压缩流动的有效方法。作者本学期学习了研究生计算流体课程，为了熟悉计算流体的基本方法，选择使用涡量-流函数法计算不可压缩方腔驱动流问题，并且对于不同雷诺数下的解进行比较和分析，得出一些结论。本文接下来的内容安排为：第 2 节提出不可压缩方腔驱动流问题，并分析该问题怎样使用涡量-流函数方法建立差分格式、选择边界条件。第 3 节介绍程序的结构。第 4 节对于不同雷诺数下的计算结果进行分析，并且与U.GHIA 等人【1】的经典结论进行对比，评述本1 1/1212文所采用的计算方法。第五节给出结论。2 2 问题的提出和分析问题的提出和分析2.12.1 经典方腔驱动流问题经典方腔驱动流问题考虑如下图所示的长度为 1 的正方形腔体，腔体上有一平板以速度 U=1 运动，其它三边为固壁条件。图 1.方腔驱动流示意图顶盖方腔驱动流问题是个很经典的问题，常常用于验证不可压缩流动数值方法的正确性。U.GHIA 等人于 1982 年发表的一篇文献（见文献【1】）计算了Re 从 100 到104的流动结果，其结果得到广泛的认同。2.22.2 涡量涡量-流函数方法简介流函数方法简介涡量-流函数法的基本思想是引入涡量和流函数：引入涡量，可以消去方程中的压力项，而引入流函数，可以使连续方程自然满足。下面对该方法进行简单推导：考虑二维问题，将式(1.2)写成分量形式：uv 0 (1.3)ttuuuP1 2u2u uv (1.4)txyxRex2y2vvvP1 2v2v uv (1.5)txyyRex2y2式(1.4)对y求偏导数减去式(1.5)对x求偏导数，考虑到=uv，推导出涡量满足yx2 2/1212122uv的方程为22(1.6)txyRexy然后引入流函数，定义为v,u(1.7)xy可见，连续性方程(1.3)自然成立。与的关系为222(1.8)2xy式(1.6)(1.8)构成了一个封闭的方程组，由(1.6)计算出涡量，再由(1.8)式计算出流函数，利用(1.7)式计算出速度。这个方程组的特点是求解速度的时候完全不用考虑压力项。若还需要求解压力场，则可以把式(1.4)对x求偏导数，式式(1.5)对y求偏导数，二者求和后整理得到关于压力的 Poisson方程22u vvu2P2(1.9)y xyx以上推导出的涡量-流函数法在计算二维问题时很成功，但是三维流动的流函数没有直观的物理意义，无法像二维流动一样直接定义，需要引入多个流函数，相应解多个 Poisson方程，计算量很大，并不实用。对于本文的二维问题，该方法就简单易行。2.32.3建立差分格式建立差分格式2.3.12.3.1划分网格划分网格方腔驱动流的流动区域很简单，均匀划分为正方形的结构网格即可，存储网格时，x 方向使用标号 i表示，y 方向使用标号 j表示，x 和 y 方向的最大网格点标号分别为 M 和 N。对于 Re 小于等于 1000 的情况，使用100*100 网格，Re 大于 1000 后的情况，使用256*256网格。计算域如图 2 所示：图 2.100*100的均分网格3 3/12122.3.22.3.2 建立差分方程建立差分方程由于本题关注的是方腔内部的流动状态，对于压力分布不关心，因此不用建立压力的差分方程。涡量的对流扩散方程(1.6)使用 FTCS 格式离散得到：1nin,ji,jt1(Reuni,jin1,jin1,j2xni,j2ni1,jvni,jin,j1in,j12yni,j2ni,j1ni1,j2(x)ni,j12(y)(1.10)该差分格式时间方向为 1 阶精度，空间方向为 2 阶精度。在(1.10)中，速度分量取的是n 时刻的值，已经对方程进行了线性化处理。流函数的 Poisson 方程中，二阶导数都用中心差分离散：11n111n1inin,in,in,1,j2ji1,jj12ji,j1x2y21in,j(1.11)这种中心差分可达到二阶精度。2.3.32.3.3 设定边界条件设定边界条件（1）速度和流函数的边界条件由于沿着壁面是一条流线，所以流函数在边界是常值，可以取为0；速度在边界满足无滑移条件。上边界（i 0 M,j N）：ui,NU 1,vi,N 0；i,N 0下边界（i 0 M,j 0）：ui,0 0,vi,0 0；i,0 0左边界（j 1 N 1,i 0）：u0,j 0,v0,j 0；0,j 0右边界（j 1 N 1,i M）：uM,j 0,vM,j 0；M,j 0（2）涡量的边界条件u2uvv2；在左右边，在上下边界，0，所以=根据涡量的定义=yyyxxv2u 0，所以=2。界，；xxy0,j=左边界（j 1 N 1,i 0）：1,j20,j1,jx2,这里引入了虚拟网格点（-1,j），4 4/12120,j 021,j=注意到，所以。1,j1,j0,j2x 0u0,j2x同理可得，右边界（j 1 N 1,i M）：M,j=2M 1,jx2下边界（i 0 M,j 0）：i,0=2i,1y2i,N 0上 边 界（i 0 M,j N）：注 意 到，所 以i,N1i,N1u1i,N2yi,N=2i,N12yy2下面考察构造的这种涡量边界条件的精度：比如对于下边界，流函数的Taylor展开为i,1i,0yyy222!y2(i,0)(i,0)y333!y3(i,0)(i,0)O(y4)(i,0)i,1i,0y2则2yyy 2!y222y 3!y333O(y4)(i,0)(i,0)i,12i,0i,1O(y4)i,12i,0i,1y2y2O(y2)对于其他边界的精度推导是类似的。所以构造的这种涡量边界条件很有优势不仅形式简单，还能有 2 阶精度。3 3 编程计算编程计算3.13.1 程序的结构程序的结构程序流程图如下图所示：5 5/1212主程序开始Cfd_driven.f90.定义变量调用 compute_node划分网格，计算节点坐标设定速度和流函数的边界条件设定涡量、流函数和速度的初始值根据流函数的值更新涡量边界条件，调用函数 compute_omiga 计算涡量场根据计算出的涡量场，调用函数compute_fai 计算流函数场输出结果根据计算出的流函数场，调用函数 compute_speed 计算速度场程序结束是是 否 收 敛？（根据相邻两次速度场之 1-范数）否图 3.流程图6 6/12123.23.2 程序说明程序说明时间步长选择 0.001，足以满足稳定性条件。1、对于涡量场的计算：根据差分方程1nin,ji,jtuni,jin1,jin1,j2xvni,jin,j1in,j12ynnnnnn1i1,j2i,ji1,ji,j12i,ji,j1()22Re(x)(y)该式中，速度分量取的是 n 时刻的值，已经对方程进行了线性化处理。所以直接使用显1示法，一步就可以将in,j求出。1nin,ji,j1in1,j2in,jin1,jin,j12in,jin,j1nin1,jin1,jin,j1in,j1nt()ui,jvi,j22Re(x)(y)2x2y（1.12）2、对于流函数场的计算：流函数满足泊松方程，差分形式为：11n111n1inin,in,in,1,j2ji1,jj12ji,j1x2y21in,j求解时要使用 JACOBI 迭代，由于x y，可以写成1,(k1)in,j1n1,(k)1,(k)n1,(k)n1,(k)2n1i1,jinhi,j1,ji,j1i,j14n+1,(k+1)n+1,(k+1)其中（k+1）表示第 k+1 次迭代，当 将此时的 n+1,(kn+1,(k+1)+1)n+1,(k)n+1,(k)时，认为精度达到，退出迭代。1n+1n+1作为 n+1 时刻的流函数场。3、速度场的求解由于 v,u，所以xy,v u i,j1i,j12yi1,ji1,j2x7 7/12124 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 接着增大，在现有二次涡的地方会出现更多的二次涡。(a)Re=100(b)Re=400(c)Re=10008 8/1212(d)Re=3200(e)Re=5000(f)Re=7500图 4.不同 Re 时的流线图(d)Re=1e4（右图出自文献 1）图 5 显示了不同 Re 数中，在通过方腔几何中心的竖直线上速度u 的分布情况，从图 5中看出，随着Re 增大，边界层越来越薄，Re5000 后，边界层变薄的速率就很小了。高Re数时，方腔中部的速度分布几乎是线性的。Re3200 时，在靠近y=1 处，u 的分布会出现一个小凸起，这个现象也被之前的文献提到过。9 9/1212图 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 中的优化网格。1010/1212图 6.不同 Re 时竖直中心线上 u 的分布对比高 Re 时，计算时间较长，如下表所示。ReRe0000CPUCPU时时间（间（s s）01111/12121 14004001000100032003200500050000 075075000000011050910109166208810940155 5 结论结论本文通过涡量-流函数法对不同雷诺数下的方腔驱动流进行了模拟，得到结论如下：对于方腔驱动流：1、随着 Re 增加，二次涡逐个出现。2、随着 Re 增加，主涡的中心逐渐移动到方腔的几何中心，Re3200 后，主涡的位置几乎不变。3、随着 Re 增加，边界层越来越薄。高Re 数时，方腔中部的速度分布几乎是线性的，在靠近 y=1 处，u 的分布会出现一个小凸起。对于本文使用的算法1、Re3200 时，计算结果精良，耗时短。2、由于网格是均分，高Re 时计算结果在边界处有较大误差。3、可以改进计算流函数泊松方程时的方法，比如使用共轭梯度法，使得收敛加快，节省计算时间。参考文献：1GHIA U,GHIA K N,SHIN C T.High-Re solutions for incompressible flowusing the Navier-Stokes Equations and a multigrid method J.J Comp utPhys,1982,48:3872411.后记：本学期学习了李嵩老师开设的计算流体课程，受益颇多。对于计算流体的基本方法和思路有所了解，于是就编写了简单的程序作为大作业。由于第一次大作业的选题过于难，没有做出，跟李老师沟通后，老师非常理解，宽限我一些时间，在此衷心感激老师的理解和宽容，就老师一学期的谆谆教诲表示感谢。1212/1212
