论文部分内容阅读
提要: 谱单元方法是一种高效的高精度计算流体动力学数值计算方法，目前被广泛运用于空气动力学的大规模模拟中。本文详细介绍了该数值计算方法好核心思想和编程思路，并实现了其程序开发。最后以单圆柱绕流问题为例验证了其准确性和高效性。模拟结果表明谱单元方法是在科学研究和工程计算中极具发展和应用前景的数值计算工具。
关键词： 谱单元、有限元、计算流体动力学、圆柱绕流
中图分类号： 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，即∂ux/∂x=0.0，∂uy/∂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=200
图4 不同雷诺数下的瞬态流线图
Fig 4 snapshots of instantaneous streamlines with different Re
从图4可以看出，对低雷诺数均匀流中圆柱绕流的数值模拟的结果与目前众多研究人员得到的结论基本上是一致的。在雷诺数Re=40时，在圆柱尾流中紧贴圆柱背后形成一对稳定的对称附着涡，没有出现漩涡脱落；随着雷诺数继续增大，稳定的对称附着涡破坏，在雷诺数Re=60附近，圆柱尾流开始出现漩涡脱落；再进一步增大雷诺数可以发现圆柱尾流中出现了成两排周期性摆动和交错的漩涡，即卡门涡街，同时还可以发现圆柱尾流初始漩涡随着雷诺数的增大逐渐向圆柱后端点靠近，流动变得更加复杂。
2.3 圆柱表面受力特性随雷诺数的变化
图5 平均阻力系数随雷诺数的变化
Fig 5 Variation of mean drag coefficients with Re
图5为平均阻力系数随雷诺数的变化，平均阻力系数随着雷诺数的增加而减小，在Re=140附近取得极小值后，随着雷诺数的增大而缓慢增大。
(a) (b)
图6 (a)阻力系数均方根 (b)升力系数均方根随雷诺数的变化
Fig 6 Variation of (a) the RMS. drag coefficients and (b) the RMS. lift coefficients with Re
图6(a)和图6(b)分别为单圆柱阻力系数均方根和升力系数均方根随雷诺数的变化，由图可知阻力系数均方根和升力系数均方根均随雷诺数的增大而逐渐增大。随着雷诺数的增大，从圆柱上脱落的涡强度增强，涡交替从圆柱上下表面脱落产生的上下表面压力脉动越来越强，因为升力系数的脉动也越来越强。
三、 结论
本文介绍了谱单元方法的插值函数选取和谱单元方法的离散，并应用于不同低雷诺数单圆柱绕流的数值模拟，得出以下结论：
（1）在雷诺数区间[40,200]，当Re=40时，圆柱后面有一对稳定对称的附着涡；在Re=60附近，圆柱后面开始出现漩涡脱落；再进一步增大雷诺数可以发现圆柱尾流中出现了成两排周期性摆动和交错的漩涡，同时还发现圆柱尾流初始漩涡随着雷诺数的增大逐渐向圆柱后端点靠近，流动变得更加复杂；Re=200的时候，圆柱后的尾流由层流向湍流转变。
（2）平均阻力系数随着雷诺数的增大而减小，在Re=140附近取得极小值后，随着雷诺数的增大而缓慢增大。阻力系數均方根和升力系数均方根均随雷诺数的增大而逐渐增大，圆柱表面压力脉动逐渐增强。
以上结果与既有文献的实验与数值模拟结果吻合，从而验证了该算法的可靠性，谱单元方法可以在较粗糙的网格下得到高精度的计算结果，表明本方法是一种有应用前景的方法。
参考文献：
[1] Gottlieb D, Orgszag S A. Numerical Analysis of Spectral Method: Theory and Application. CBMs-NFS Monograph No 26 [M]. Society for Industry and Applied Mathematics Philadelphia, 1977.
[2] Gottlieb D, Hussaini M Y, Orgszag S A. Theory and Application of Spectral Method [A]. In: Voigt R G. Gottlieb D. Hussaini M Y. Spectral Method for Partial Differential Equation [C]. SIAM Philadelphia, 1984: 1-54.
[3] Canuto C, Quarteroni A, Hussaini M Y, Zang T A. Spectral methods in fluid dynamics [M]. Springer series in computational physics, spriger-verlag New York Inc. 1988.
[4] Patera A T. A spectral element method for fluid dynamics: laminar flow in a channel expansion [J]. J Comput Phys, 1984, 154: 468-488.
[5] Korczak K E, Patera A T. An isoparametric spectral element method for solution of the Navier-Stokes equation in complex geometry [J]. Journal of Computer Physics, 1986, 62: 361-382.
[6] 韩庆书, 王唯. 求解二维Navier-Stokes方程的谱元法 [J]. 北京联合大学学报, 1998, 12(S1-S2): 123-131.
[7] Gear C W. Numerical initial value in ordinary differential equations [M]. Prentice-Hall Englewwod Cliffs, NJ, 1973.
[8] 陈雪江, 秦国良, 徐忠. 谱元法和高阶时间分裂法求解方腔顶盖驱动流 [J].计算力学学报, 2002, 19(3): 281-285.
关键词： 谱单元、有限元、计算流体动力学、圆柱绕流
中图分类号： 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，即∂ux/∂x=0.0，∂uy/∂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=200
图4 不同雷诺数下的瞬态流线图
Fig 4 snapshots of instantaneous streamlines with different Re
从图4可以看出，对低雷诺数均匀流中圆柱绕流的数值模拟的结果与目前众多研究人员得到的结论基本上是一致的。在雷诺数Re=40时，在圆柱尾流中紧贴圆柱背后形成一对稳定的对称附着涡，没有出现漩涡脱落；随着雷诺数继续增大，稳定的对称附着涡破坏，在雷诺数Re=60附近，圆柱尾流开始出现漩涡脱落；再进一步增大雷诺数可以发现圆柱尾流中出现了成两排周期性摆动和交错的漩涡，即卡门涡街，同时还可以发现圆柱尾流初始漩涡随着雷诺数的增大逐渐向圆柱后端点靠近，流动变得更加复杂。
2.3 圆柱表面受力特性随雷诺数的变化
图5 平均阻力系数随雷诺数的变化
Fig 5 Variation of mean drag coefficients with Re
图5为平均阻力系数随雷诺数的变化，平均阻力系数随着雷诺数的增加而减小，在Re=140附近取得极小值后，随着雷诺数的增大而缓慢增大。
(a) (b)
图6 (a)阻力系数均方根 (b)升力系数均方根随雷诺数的变化
Fig 6 Variation of (a) the RMS. drag coefficients and (b) the RMS. lift coefficients with Re
图6(a)和图6(b)分别为单圆柱阻力系数均方根和升力系数均方根随雷诺数的变化，由图可知阻力系数均方根和升力系数均方根均随雷诺数的增大而逐渐增大。随着雷诺数的增大，从圆柱上脱落的涡强度增强，涡交替从圆柱上下表面脱落产生的上下表面压力脉动越来越强，因为升力系数的脉动也越来越强。
三、 结论
本文介绍了谱单元方法的插值函数选取和谱单元方法的离散，并应用于不同低雷诺数单圆柱绕流的数值模拟，得出以下结论：
（1）在雷诺数区间[40,200]，当Re=40时，圆柱后面有一对稳定对称的附着涡；在Re=60附近，圆柱后面开始出现漩涡脱落；再进一步增大雷诺数可以发现圆柱尾流中出现了成两排周期性摆动和交错的漩涡，同时还发现圆柱尾流初始漩涡随着雷诺数的增大逐渐向圆柱后端点靠近，流动变得更加复杂；Re=200的时候，圆柱后的尾流由层流向湍流转变。
（2）平均阻力系数随着雷诺数的增大而减小，在Re=140附近取得极小值后，随着雷诺数的增大而缓慢增大。阻力系數均方根和升力系数均方根均随雷诺数的增大而逐渐增大，圆柱表面压力脉动逐渐增强。
以上结果与既有文献的实验与数值模拟结果吻合，从而验证了该算法的可靠性，谱单元方法可以在较粗糙的网格下得到高精度的计算结果，表明本方法是一种有应用前景的方法。
参考文献：
[1] Gottlieb D, Orgszag S A. Numerical Analysis of Spectral Method: Theory and Application. CBMs-NFS Monograph No 26 [M]. Society for Industry and Applied Mathematics Philadelphia, 1977.
[2] Gottlieb D, Hussaini M Y, Orgszag S A. Theory and Application of Spectral Method [A]. In: Voigt R G. Gottlieb D. Hussaini M Y. Spectral Method for Partial Differential Equation [C]. SIAM Philadelphia, 1984: 1-54.
[3] Canuto C, Quarteroni A, Hussaini M Y, Zang T A. Spectral methods in fluid dynamics [M]. Springer series in computational physics, spriger-verlag New York Inc. 1988.
[4] Patera A T. A spectral element method for fluid dynamics: laminar flow in a channel expansion [J]. J Comput Phys, 1984, 154: 468-488.
[5] Korczak K E, Patera A T. An isoparametric spectral element method for solution of the Navier-Stokes equation in complex geometry [J]. Journal of Computer Physics, 1986, 62: 361-382.
[6] 韩庆书, 王唯. 求解二维Navier-Stokes方程的谱元法 [J]. 北京联合大学学报, 1998, 12(S1-S2): 123-131.
[7] Gear C W. Numerical initial value in ordinary differential equations [M]. Prentice-Hall Englewwod Cliffs, NJ, 1973.
[8] 陈雪江, 秦国良, 徐忠. 谱元法和高阶时间分裂法求解方腔顶盖驱动流 [J].计算力学学报, 2002, 19(3): 281-285.