基于谱单元方法的单圆柱绕流特性分析.docVIP

基于谱单元方法的单圆柱绕流特性分析 　　提要: 谱单元方法是一种高效的高精度计算流体动力学数值计算方法，目前被广泛运用于空气动力学的大规模模拟中。本文详细介绍了该数值计算方法好核心思想和编程思路，并实现了其程序开发。最后以单圆柱绕流问题为例验证了其准确性和高效性。模拟结果表明谱单元方法是在科学研究和工程计算中极具发展和应用前景的数值计算工具。 　　关键词： 谱单元、有限元、计算流体动力学、圆柱绕流 　　中图分类号： O313 文献标识码： A 文章编号： 　　自从1977年Gottlieb和Orszag[1]系统地从数学方面对谱方法进行了理论的阐述，它与有限差分法及有限元法一起构成了求解偏微分方程的三大方法，被广泛地应用于更多的领域。随着谱方法在各领域的应用和发展，谱方法在理论研究上日趋完善，它开辟了谱方法应用函数分析技术处理复杂问题的道路。1984年，Gottlieb和Hussaini开始将谱方法向计算流体动力学方面推广[2,3]。到了80年代初期，Patera才结合谱方法的精度和有限元的思想提出所谓的谱单元方法[4]，谱单元方法具有谱方法的高精度和收敛特性，并且还可以像有限元法一样具有很好的几何区域的适应性[5]。 　　本文研究了谱单元方法插值函数的选取和谱单元的离散过程，给出了离散方程的一般形式，并采用时间分裂格式的谱单元法求解Navier-Stokes方程，以不同雷诺数下单圆柱绕流的数值模拟作为基本算例，验证了谱单元法的高精度和计算效率，计算表明结果令人满意。 　　 　　一、 谱单元离散格式 　　 　　 　　二、 单圆柱绕流计算分析 　　在研究圆柱流场时常用的几个无量纲化系数：CD（阻力系数）,CL（升力系数）和 St（斯托罗哈数）定义如下： 　　 （12） 　　其中，FD为阻力，与来流方向一致，主要由流体绕圆柱柱表面摩擦阻力以及圆柱前后压力差造成；FL为升力，与来流方向垂直，主要由涡交替从圆柱上下表面脱落产生上下表面压力脉动造成；St为涡脱落频率，D为圆柱直径。 　　2.1 计算域和网格划分 　　考虑直径为D的圆柱受到未经扰动的均匀来流作用，基于圆柱直径和来流流速的雷诺数取Re=200。所选计算域50D×40D，圆柱位于坐标系原点（0,0）。入口边界和出口边界分别位于圆柱中心上游20D和下游30D处，流域顶部和底部离圆柱中心20D。相应的边界条件如下：进口处自由来流速度为绕流问题特征速度，即ux=U∞，uy=0.0；上下边界条件与进口边界条件相同；出口边界处纵向和横向速度梯度均为0.0，即#8706;ux/#8706;x=0.0，#8706;uy/#8706;x=0.0；圆柱表面处为不可滑移边界条件，即ux=0.0，uy=0.0。计算域和边界条件如图1所示。 　　 　　 　　图1 计算域和边界条件示意图 　　Fig 1 Schematic diagram of the computational domain and boundary conditions 　　 　　计算域网格划分采用了四边形非结构化谱单元网格，总共划分了354个单元，如图2(a)所示。在靠近圆柱壁面的地方进行了几层非常细的网格加密，离圆柱壁面最近的一层网格厚度为0.1D，如图2(b)所示。同时，在圆柱尾流区域也进行了加密处理。 　　 　　图2 (a)谱单元网格划分示意图 (b)圆柱附近网格加密示意图 　　Fig 2 (a) spectral element mesh, 354 elements (b) zoomed-in view of the mesh around the cylinder 　　 　　为了验证插值函数的阶数对计算结果的影响，对单圆柱绕流进行了基于三种不同阶数的插值函数的数值模拟。在算例1中，谱函数插值采用了N=5阶GLL二维拉格朗日形函数；在算例2中，N=7；在算例3中，N=9。计算时间步长为Δt=0.005。如图3所示为所得阻力系数和升力系数时程曲线。 　　 　　 　　 (a)(b) 　　图3 单圆柱绕流阻力系数(a)和升力系数(b)时程曲线 　　Fig 3 Time histories of drag and lift coefficients for a cylinder 　　 　　从图中可以看到，雷诺数为200时单圆柱绕流的阻力系数与升力系数均呈周期性正弦变化，而且升力系数变化周期是阻力系数变化周期的两倍，这是由于漩涡交替从圆柱上下表面脱落。还可以看到，N=7和N=9的曲线几乎吻合在一起，极为相似。 　　 　　2.2 单圆柱绕流的流态 　　 　　 　　(a) Re=40 (b) Re=60 　　 　　(a) Re=120 (b) Re