? 带横隔板圆柱绕流特性数值模拟 带横隔板圆柱绕流特性数值模拟
水庆象， 王大国，王 伟
(西南科技大学环境与资源学院，四川绵阳 621010)
摘要：为揭示尾迹区添加横隔板对圆柱绕流流场特性影响，把多步格式引入到特征线算子分裂有限元法中，建立了基于多步格式的特征线算子分裂有限元法：在每个时间步内将Navier-Stokes方程分裂成对流项和扩散项，对流项时间离散采用多步格式，在每一子时间步内沿特征线展开并显式求解。方腔流数值模拟结果表明该算法既可降低对整体时间步长的要求又可提高计算精度。对比有无横隔板圆柱绕流流场和圆柱表面压力变化表明，横隔板可以有效地抑制绕流尾迹区涡旋脱落，提高圆柱背流面压力，减少圆柱上下表面的压力差。
关键词：多步格式；有限元；方腔流； 圆柱绕流；横隔板
钝体绕流常见于水利工程、土木工程等领域，其中圆柱绕流应用较为广泛。当流体绕过这些圆柱形结构物时会在其表面发生边界分离，在一定雷诺数下圆柱后方出现流动分离和涡脱[1]。随着涡脱的产生，旋涡会对结构物产生顺着来流方向(阻力)和垂直来流方向(升力)周期变化的作用力，这不仅会诱发结构物的疲劳损伤，甚至还会引起共振，影响结构物的可靠性和安全性，缩短结构物寿命。因此，如何通过合适的措施减小结构物的阻力、抑制其流动分离进而减小周期性的脉动升力已成为工程上亟待解决的问题。圆柱尾迹区旋涡形成和脱落的机理复杂，Gerrard[2]指出上下剪切层间相互作用是旋涡形成和脱落的关键。据此，学者们提出在圆柱后添加横隔板来限制旋涡脱落，改变升阻力系数。Kwon和Choi[3]对Re分别为100、160的带不同长度横隔板圆柱绕流展开研究，给出了平均阻力系数和斯特劳哈数的变化趋势，两种Re下平均阻力系数随着横隔板长度的增加而降低，而斯特劳哈数则有着不同的变化趋势；但没有对升力系数的变化展开研究。段志强[4]采用二维层流模型研究Re=200的带横隔板圆柱绕流，虽定性地给出了升力系数随横隔板长度变化的趋势，但发现横隔板长度与圆柱直径相等时升力系数幅值会减小为0；而Kawai[5]的研究表明，当横隔板长度与圆柱直径相等时升力系数幅值会减小但不会为0。
研究带横隔板圆柱绕流场特性的方法主要有实验方法和数值模拟方法。数值模拟方法研究带横隔板圆柱绕流问题的关键是求解Navier-Stokes(N-S)方程。在求解具有复杂几何边界流动问题时，有限元法更具优势。但当雷诺数较大使得对流占优时，N-S方程呈现高度的非线性特性，采用标准Galerkin有限元法求解将引起数值解的振荡。针对对流项占优时导致数值解振荡的问题，已有多种有别于标准Galerkin有限元法的新算法问世，如Petrov-Galerkin法[6]、Taylor-Galerkin法[7]、Galerkin最小二乘法[8]、特征线Galerkin法[9]等。其中，特征线Galerkin法的核心是将N-S方程的动量方程沿特征线进行时间离散，以理性形式引入稳定项，具有较明确的物理意义。1995年，Zienkiewicz和Codina[10]将Taylor展开法引入到特征线Galerkin法中并将其与分离算法相结合，提出了基于特征线的分离算法(CBS法)，该法可直接由N-S方程推导出合理的平衡耗散项，可获得稳定的数值解。Wang等[11]结合CBS法和算子分裂法[12]的优点，提出特征线算子分裂有限元法(CBOS有限元法)求解N-S方程。该方法在处理扩对流项时采用显式格式，显式格式具有单个时间步的计算量小、程序易于实现等优点，但其不足为在时域离散化计算中是条件稳定的[13]，即必须采用较小的时间步长才能使计算稳定。因此，对流项的显式求解使得CBOS有限元法受时间步长的限制，影响计算效率。
本文把多步格式引入到CBOS有限元法中，提出基于多步格式的特征线算子分裂有限元法。通过方腔流验证其优越性后，将该方法用于带与圆柱直径等长的横隔板圆柱绕流，以揭示尾迹区放置横隔板对圆柱绕流流场参数的影响，并结合流场和圆柱表面压力的变化阐明添加横隔板对圆柱绕流流动参数的影响机理。
1 数值模型
1.1 控制方程
二维非定常不可压缩流体基本控制方程的量纲一形式表示为
(1)
式中： Ω和
〗分别为空间域和时间域； i,j=1,2；(ui,uj)=(u,v)，u和v分别为水平方向速度和垂直方向速度；p为压力；t为时间；
为雷诺数(其中μ为动力黏度，ρ为流体密度，U为特征速度，L为特征长度)；fi为外力；(x1,x2)=(x,y)。
1.2 CBOS有限元法
在每一个时间步内，采用算子分裂法将N-S方程分为扩散项和对流项：
(2)
(3)
式中：
为扩散项(式(2))在n+1时刻的解，同时为对流项(式(3))在n+1时刻的初值；
为对流项在n+1时刻的解，同时也是N-S方程式(1)在n+1时刻的解。
式(2)和式(3)应用CBOS有限元法[11]，可得如下时间离散形式：
(4)
(5)
式(5)为对流项式(3)沿特征线显式时间离散所得。
1.3 基于多步格式的CBOS有限元法
采用多步格式处理对流项式(3)，即在每一个整体时间步Δt内将对流项分为h步进行求解，则求解对流项的时间步长为
。根据式(5)，对流项在第l子时间步内的求解公式为
(6)
式中： l依次取1,2,…,h。若l=1，则
；若l=h，则
。
对式(4)和式(6)采用标准Galerkin加权余量有限元法进行空间离散。基于多步格式的CBOS有限元法的求解过程为：
(1)以n时刻的速度场
和压力场pn作为初值，求解式(2)，得到n+1时刻速度场的过渡值
和压力场pn+1；
(2)以
为初值，由式(3)求解速度场
。在求解式(3)时，采用多步格式，在每一子时间步内沿特征线展开得式(6)并显式求解；
(3)转到下一时刻，重复步骤(1)和步骤(2)。
2 模型验证
在式(6)中，若取h=1，即对流项不采用多步格式，此时为CBOS有限元法。为比较本文方法与CBOS有限元法(h=1时)的不同数值表现，分别采用这两种算法对方腔流进行数值模拟。方腔流几何模型及边界条件如图1(a)所示。取L=1 m，U=1 m/s，ρ=1.0 kg/m3，调整动力黏度得到不同的雷诺数。左下角为坐标原点，在该点相对压力为0。整个计算域划分为30×30个九节点四边形单元，共有3 721个节点，具体网格剖分如图1(b)。
图1 方腔流
Fig.1Driven cavity flow
取动力黏度μ=0.001 Pa·s， 即Re=1 000，分别取Δt=0.1、0.05、0.01和h=1、5、10。当h=1时即采用CBOS有限元法进行数值模拟，Δt分别取0.1和0.05时计算值发散，需要取Δt=0.01时才能收敛；当h分别取5和10时，分别采用Δt为0.1、0.05和0.01均可收敛。
图2给出了Δt=0.01时不同对流子时间步下速度沿中线的分布。从图2可以看出，Δt=0.01时无论h如何取值计算均收敛，但h=1时计算精度较低，而h增大后精度得到提高；当h=10时，本文计算结果与Ghia等[14]结果吻合较好。由以上分析可知CBOS有限元法由于对流项的显式求解使得该算法受到时间步长的束缚，而当对流项采用多步格式既可降低对整体时间步长的要求又可提高计算精度。
图2 Δt=0.01时不同对流子时间步下速度沿中线的分布
Fig.2Velocity along lines through geometric center at the different h when Δt=0.01
3 带横隔板圆柱绕流特性分析
3.1 计算区域和边界条件
图3给出了带横隔板圆柱结构尺度示意图及其计算网格剖分图，LS为横隔板长度，R为圆柱半径。计算中取LS=2R(带与圆柱直径等长的横隔板)，横隔板厚度为
。取L=1 m，U=1 m/s， ρ=1.0 kg/m3，调整动力黏度得到不同的雷诺数。入口处指定沿水平方向的均匀来流u=1，垂向速度v=0；侧壁、圆柱表面和横隔板表面为固壁边界条件；出口为自由出流边界，在该位置相对压力为0。对流子时间步取h=10。
图3 带横隔板圆柱绕流结构及其计算网格剖分示意
Fig.3Size diagram and sketch of computational grid for a circular cylinder with splitter plate
3.2 结果分析
3.2.1 流场参数特性分析
表1列出了本文计算的圆柱绕流(LS =0)和带横隔板圆柱绕流平均阻力系数
、升力系数CL和斯特劳哈数St。从表1可知，在圆柱背流侧增加横隔板，当Re=100时，平均阻力系数
由1.394减小到1.223，减少了约12.3%；升力系数由0.332减少为0.076，减小了约77.1%；St由0.167减小为0.138，减少了约17.37%；在Hwang等[19]的研究中分别减少了约12.7%、76.8%和17.96%。当Re=200时，平均阻力
由1.339减小到1.0624，减少了约20.6%；圆柱所受升力幅值由0.614减少为0.163，减小了约73.5%；St由0.194减小为0.155，减少了约20.1%。计算结果表明二维层流状态下在圆柱尾迹区内放置与圆柱直径等长的横隔板使得平均阻力系数、升力系数幅值和斯特劳哈数都会减小，但都不会减为0。而段志强[4]采用二维层流模型研究Re=200的带横隔板圆柱绕流，发现横隔板长度与圆柱直径相等时升力系数幅值会减小为0。
表1 Re=100和Re=200时本文计算结果
Table 1Flow parameters for flow field at Re=100 and Re=200
参数不带横隔板(LS=0)带横隔板Re=100Re=200Re=100Re=200本文结果Han等[15]Xu和Wang[16]本文结果Ding等[17]Harichandan和Roy[18]本文结果本文结果CD1.3941.3791.4231.3391.3481.3301.2321.062CL±0.332±0.334±0.340±0.614±0.659±0.600±0.076±0.163St0.1670.1680.1710.1940.1960.1920.1380.155
3.2.2 绕流场结果特性分析
图4和图5分别给出了Re=200时圆柱绕流和带横隔板圆柱绕流半个周期内的涡量图。图4和图5对比可见：由于横隔板的存在，圆柱壁面附近的自由剪切层不能相互作用，剪切层只能沿着流动方向向下游输送，直至到达横隔板末端后，上下剪切层才开始相互作用。圆柱尾迹区的旋涡变长，相应涡脱周期变长，St变小，说明在圆柱尾部添加横隔板能有效地抑制圆柱背流面旋涡脱落。
图4 Re=200时圆柱绕流半个周期内涡量
Fig.4Vortices distribution for a circular cylinder at Re=200 during half cycle
图5 Re=200时带横隔板圆柱绕流半个周期内涡量
Fig.5Vortices distribution for a circular cylinder with splitter plate at Re=200 during half cycle
图6给出了Re=200时的圆柱绕流在大致一个涡旋脱落周期内的压力和流线的演化过程。图7给出了Re=200时带横隔板圆柱绕流一个周期内压力和流线演化。图6和图7对比可知，在两种情况下圆柱迎流面压力基本不变，而带横隔板圆柱背流面压力较不带横隔板圆柱背流面压力有着显著的增加，导致带横隔板圆柱所受平均阻力系数降低。
图6 Re=200时一个周期内压力和流线演化
Fig.6Development of pressure and streamline during a cycle at Re=200
图7 Re=200时带横隔板圆柱绕流一个周期内压力和流线演化
Fig.7Development of pressure and streamline for a circular cylinder with splitter plate during a cycle at Re=200
由图4和图5对比可见，带横隔板圆柱绕流流场尾迹区内涡量分布的对称性较不带横隔板圆柱绕流情况有了比较明显的改善，流场对称性的改善可改变圆柱上的升力系数幅值。为揭示升力幅值减小且不减为零的原因，图8给出了Re=200正向升力系数达到最大时圆柱表面压力系数Cp分布。由图8可以看出，当正向升力系数达到最大时，带横隔板圆柱背流面上部(90°
(实线：不带横隔板圆柱绕流计算结果；虚线：带横隔板圆柱绕流计算结果)
图8 Re=200正向升力系数达到最大时圆柱表面压力系数分布
Fig.8Distribution of surface pressure coefficient when the maximum of the positive lift coefficient is reached at Re=200
4 结 论
(1)CBOS有限元法中显式求解对流项在时域离散化计算中是条件稳定的，把多步格式引入到特征线算子分裂有限元法中，建立了基于多步格式的特征线算子分裂有限元法，该算法中对流项时间离散采用多步格式，在每一子时间步内沿特征线展开并显式求解。方腔流数值模拟结果证实该算法既可降低对整体时间步长的要求又可提高计算精度。
(2)带与圆柱直径等长的横隔板圆柱绕流模拟结果表明：当Re=100时，平均阻力系数
减少了12.3%；升力系数幅值减小了77.1%；St减少了17.37%；当Re=200时，3个流场参数分别减少20.6%、73.5%和20.1%。
(3)通过对比Re=200时有无横隔板绕流流场和圆柱表面压力变化表明：圆柱尾迹区放置横隔板后，剪切层只能沿着流动方向向下游输送，直至到达横隔板末端后，上下剪切层才开始相互作用，旋涡变长，旋涡脱落受到抑制，涡脱周期变长，St变小；带横隔板圆柱背流面压力较不带横隔板圆柱背流面压力有着显著的增加而迎流面压力无明显变化，导致带横隔板圆柱所受平均阻力系数降低；横隔板的加入使得圆柱上下表面压力系数差大幅度减小但仍存在差距，所以带横隔板圆柱所受升力系数幅值也大幅度降低但不会减少至0。
参考文献：
[1]吴剑, 齐鄂荣, 李炜, 等. 应用PIV系统研究横流中近壁水平圆柱绕流旋涡特性[J]. 水科学进展, 2005, 16(5): 628-633. (WU J, QI E R, LI W, et al. Research on vortexes of cross flow past a horizontal circular cylinder by PIV system [J]. Advances in Water Science, 2005, 16(5): 628-633. (in Chinese))
[2]GERRARD J H. The mechanics of the formation region of vortices behind bluff bodies [J]. Journal of Fluid Mechanics, 1966, 25(2): 401-413.
[3]KWON K, CHOI H. Control of laminar vortex shedding behind a circular cylinder using splitter plates[J]. Physics of Fluids, 1996, 8(2): 479-486.
[4]段志强. 低雷诺数下尾部隔板影响圆柱绕流的数值研究[D]. 重庆: 重庆大学, 2012: 16-25. (DUAN Z Q. An numerical investigation on fluid flow past a circular cylinder with splitter plate at low Reynolds number [D]. Chongqing: Chongqing University, 2012: 16-25. (in Chinese))
[5]KAWAI H. Discrete vortex simulation for flow around a circular cylinder with a splitter plate[J]. Journal of Wind Engineering and Industrial Areodynamics, 1990, 33(1/2): 153-160.
[6]SATAPRAHM C, LUADSONG A. The meshless local Petrov-Galerkin method for simulating unsteady incompressible fluid flow [J]. Journal of the Egyptian Mathematical Society, 2014, 22(3): 501-510.
[7]KUMAR B V R, KUMAR S. Convergence of three-step Taylor Galerkin finite element scheme based monotone Schwarz iterative method for singularly perturbed differential-difference equation [J]. Numerical Functional Analysis and Optimization, 2015, 36(8): 1029-1045.
[8]BARBOSA A G, DUTRA E G, DUDA F P. A sensitivity analysis on the parameter on the GLS for a second-gradient theory of incompressible flow [J]. International Journal for Numerical Methods in Fluids, 2015, 78(10): 638-655.
[9]TAN L, ZHU B S, WANG Y C, et al. Turbulent flow simulation using large eddy simulation combined with characteristic-based split scheme [J]. Computers& Fluids, 2014, 94(1): 161-172.
[10]ZIENKIEWICZ O C, CODINA R. A general algorithm for compressible and incompressible-flow: I: the split characteristic based scheme [J]. International Journal for Numerical Methods in Fluids, 1995, 20(8/9): 869-885.
[11]WANG D G, WANG H J, XIONG J H, et al. Characteristic-based operator-splitting finite element method for Navier-Stokes equations [J]. Science China (Technological Sciences), 2011, 54(8): 2157-2166.
[12]GLOWINSKI R, PIRONNEAU O. Finite element method for the Navier-Stokes equations [J]. Annual Review of Fluid Mechanics, 1992, 24(1): 167-204.
[13]韩先洪, 李锡夔. 高黏性流动有限元模拟的迭代稳定分步算法 [J]. 力学学报, 2006, 38(1): 16-24. (HAN X H, LI X K. An iterative stabilized fractional step algorithm for finite element analysis in high-viscosity fluid flows[J]. Chinese Journal of Theoretical and Applied Mechanics, 2006, 38(1): 16-24. (in Chinese))
[14]GHIA U, GHIA K N, SHIN C T. High-Re solutions for incompressible flow using the Navier-Stokes equations and a multigrid method [J]. Journal of Computational Physics, 1982, 48(3): 387-411.
[15]HAN Z L, ZHOU D, TU J H. Laminar flow patterns around three side-by-side arranged circular cylinders using semi-implicit three-step Taylor-characteristic-based-split(3-TCBS) algorithm [J]. Engineering Applications of Computational Fluid Mechanics, 2013, 7(1): 1-12.
[16]XU S, WANG Z J. An immersed interface method for simulating the interaction of a fluid with moving boundaries [J]. Journal of Computational Physics, 2006, 216(2): 454-493.
[17]DING H, SHU C, YEO K S, et al. Numerical simulation of flows around two circular cylinders by mesh-free least square-based finite difference methods [J]. International Journal for Numerical Methods in Fluids, 2007, 53(2): 305-332.
[18]HARICHANDAN A B, ROY A. Numerical investigation of low Reynolds number flow past two and three circular cylinders using unstructured grid CFR scheme [J]. International Journal of Heat and Fluid Flow, 2010, 31(2): 154-171.
[19]HWANG J Y, YANG K S, SUN S H. Reduction of flow-induced forces on a circular cylinder using a detached splitter plate [J]. Physics of Fluids, 2003, 15(8): 2433-2436.
*The study is financially supported by the National Natural Science Foundation of China (No.51349011).
Numerical simulation of the characteristics of flow pastcircular cylinder with splitter plate
SHUI Qingxiang, WANG Daguo, WANG Wei
(School of Environment and Resource, Southwest University of Science and Technology, Mianyang 621010, China)
Abstract:A multistep characteristic-based operator-splitting finite element method is proposed to determine the effects of adding a splitter plate at the back of a cylinder on the characteristics of a flow past a circular cylinder. In each time step, the N-S equations are split into diffusive and convective parts. For the convective part, temporal discretization is performed by the characteristic Galerkin method with a multistep technique. Results of the present model are in good agreement with previously published data. Thus, the present model is competitive in solving incompressible viscous flows in terms of accuracy and efficiency. Furthermore, the flow past a circular cylinder with a splitter plate was simulated by the present model. The flow fields and pressure on the cylinder surface with and without a splitter plate were compared. Findings indicated that the plate can suppress vortex shedding, increase the pressure on the leeward of the cylinder, and reduce the pressure difference between the upper and lower surfaces.
Key words:multistep-scheme; finite element method; driven cavity flow; flow past a circular cylinder; splitter plate
DOI：10.14042/j.cnki.32.1309.2016.04.013
收稿日期：2015-09-28;网络出版时间：2016-06-22
基金项目：国家自然科学基金资助项目(51349011)；四川省教育厅科研项目(16CZ0013)
作者简介：水庆象(1988—)，男，安徽蚌埠人，讲师，硕士，主要从事计算水力学方面研究。
E-mail:[email protected]
通信作者：王大国， E-mail:[email protected]
中图分类号：TV135;O35
文献标志码：A
文章编号：1001-6791(2016)04-0586-08
网络出版地址：http://www.cnki.net/kcms/detail/32.1309.P.20160622.2228.028.html
? 带横隔板圆柱绕流特性数值模拟 带横隔板圆柱绕流特性数值模拟
水庆象， 王大国，王 伟
(西南科技大学环境与资源学院，四川绵阳 621010)
摘要：为揭示尾迹区添加横隔板对圆柱绕流流场特性影响，把多步格式引入到特征线算子分裂有限元法中，建立了基于多步格式的特征线算子分裂有限元法：在每个时间步内将Navier-Stokes方程分裂成对流项和扩散项，对流项时间离散采用多步格式，在每一子时间步内沿特征线展开并显式求解。方腔流数值模拟结果表明该算法既可降低对整体时间步长的要求又可提高计算精度。对比有无横隔板圆柱绕流流场和圆柱表面压力变化表明，横隔板可以有效地抑制绕流尾迹区涡旋脱落，提高圆柱背流面压力，减少圆柱上下表面的压力差。
关键词：多步格式；有限元；方腔流； 圆柱绕流；横隔板
钝体绕流常见于水利工程、土木工程等领域，其中圆柱绕流应用较为广泛。当流体绕过这些圆柱形结构物时会在其表面发生边界分离，在一定雷诺数下圆柱后方出现流动分离和涡脱[1]。随着涡脱的产生，旋涡会对结构物产生顺着来流方向(阻力)和垂直来流方向(升力)周期变化的作用力，这不仅会诱发结构物的疲劳损伤，甚至还会引起共振，影响结构物的可靠性和安全性，缩短结构物寿命。因此，如何通过合适的措施减小结构物的阻力、抑制其流动分离进而减小周期性的脉动升力已成为工程上亟待解决的问题。圆柱尾迹区旋涡形成和脱落的机理复杂，Gerrard[2]指出上下剪切层间相互作用是旋涡形成和脱落的关键。据此，学者们提出在圆柱后添加横隔板来限制旋涡脱落，改变升阻力系数。Kwon和Choi[3]对Re分别为100、160的带不同长度横隔板圆柱绕流展开研究，给出了平均阻力系数和斯特劳哈数的变化趋势，两种Re下平均阻力系数随着横隔板长度的增加而降低，而斯特劳哈数则有着不同的变化趋势；但没有对升力系数的变化展开研究。段志强[4]采用二维层流模型研究Re=200的带横隔板圆柱绕流，虽定性地给出了升力系数随横隔板长度变化的趋势，但发现横隔板长度与圆柱直径相等时升力系数幅值会减小为0；而Kawai[5]的研究表明，当横隔板长度与圆柱直径相等时升力系数幅值会减小但不会为0。
研究带横隔板圆柱绕流场特性的方法主要有实验方法和数值模拟方法。数值模拟方法研究带横隔板圆柱绕流问题的关键是求解Navier-Stokes(N-S)方程。在求解具有复杂几何边界流动问题时，有限元法更具优势。但当雷诺数较大使得对流占优时，N-S方程呈现高度的非线性特性，采用标准Galerkin有限元法求解将引起数值解的振荡。针对对流项占优时导致数值解振荡的问题，已有多种有别于标准Galerkin有限元法的新算法问世，如Petrov-Galerkin法[6]、Taylor-Galerkin法[7]、Galerkin最小二乘法[8]、特征线Galerkin法[9]等。其中，特征线Galerkin法的核心是将N-S方程的动量方程沿特征线进行时间离散，以理性形式引入稳定项，具有较明确的物理意义。1995年，Zienkiewicz和Codina[10]将Taylor展开法引入到特征线Galerkin法中并将其与分离算法相结合，提出了基于特征线的分离算法(CBS法)，该法可直接由N-S方程推导出合理的平衡耗散项，可获得稳定的数值解。Wang等[11]结合CBS法和算子分裂法[12]的优点，提出特征线算子分裂有限元法(CBOS有限元法)求解N-S方程。该方法在处理扩对流项时采用显式格式，显式格式具有单个时间步的计算量小、程序易于实现等优点，但其不足为在时域离散化计算中是条件稳定的[13]，即必须采用较小的时间步长才能使计算稳定。因此，对流项的显式求解使得CBOS有限元法受时间步长的限制，影响计算效率。
本文把多步格式引入到CBOS有限元法中，提出基于多步格式的特征线算子分裂有限元法。通过方腔流验证其优越性后，将该方法用于带与圆柱直径等长的横隔板圆柱绕流，以揭示尾迹区放置横隔板对圆柱绕流流场参数的影响，并结合流场和圆柱表面压力的变化阐明添加横隔板对圆柱绕流流动参数的影响机理。
1 数值模型
1.1 控制方程
二维非定常不可压缩流体基本控制方程的量纲一形式表示为
(1)
式中： Ω和
〗分别为空间域和时间域； i,j=1,2；(ui,uj)=(u,v)，u和v分别为水平方向速度和垂直方向速度；p为压力；t为时间；
为雷诺数(其中μ为动力黏度，ρ为流体密度，U为特征速度，L为特征长度)；fi为外力；(x1,x2)=(x,y)。
1.2 CBOS有限元法
在每一个时间步内，采用算子分裂法将N-S方程分为扩散项和对流项：
(2)
(3)
式中：
为扩散项(式(2))在n+1时刻的解，同时为对流项(式(3))在n+1时刻的初值；
为对流项在n+1时刻的解，同时也是N-S方程式(1)在n+1时刻的解。
式(2)和式(3)应用CBOS有限元法[11]，可得如下时间离散形式：
(4)
(5)
式(5)为对流项式(3)沿特征线显式时间离散所得。
1.3 基于多步格式的CBOS有限元法
采用多步格式处理对流项式(3)，即在每一个整体时间步Δt内将对流项分为h步进行求解，则求解对流项的时间步长为
。根据式(5)，对流项在第l子时间步内的求解公式为
(6)
式中： l依次取1,2,…,h。若l=1，则
；若l=h，则
。
对式(4)和式(6)采用标准Galerkin加权余量有限元法进行空间离散。基于多步格式的CBOS有限元法的求解过程为：
(1)以n时刻的速度场
和压力场pn作为初值，求解式(2)，得到n+1时刻速度场的过渡值
和压力场pn+1；
(2)以
为初值，由式(3)求解速度场
。在求解式(3)时，采用多步格式，在每一子时间步内沿特征线展开得式(6)并显式求解；
(3)转到下一时刻，重复步骤(1)和步骤(2)。
2 模型验证
在式(6)中，若取h=1，即对流项不采用多步格式，此时为CBOS有限元法。为比较本文方法与CBOS有限元法(h=1时)的不同数值表现，分别采用这两种算法对方腔流进行数值模拟。方腔流几何模型及边界条件如图1(a)所示。取L=1 m，U=1 m/s，ρ=1.0 kg/m3，调整动力黏度得到不同的雷诺数。左下角为坐标原点，在该点相对压力为0。整个计算域划分为30×30个九节点四边形单元，共有3 721个节点，具体网格剖分如图1(b)。
图1 方腔流
Fig.1Driven cavity flow
取动力黏度μ=0.001 Pa·s， 即Re=1 000，分别取Δt=0.1、0.05、0.01和h=1、5、10。当h=1时即采用CBOS有限元法进行数值模拟，Δt分别取0.1和0.05时计算值发散，需要取Δt=0.01时才能收敛；当h分别取5和10时，分别采用Δt为0.1、0.05和0.01均可收敛。
图2给出了Δt=0.01时不同对流子时间步下速度沿中线的分布。从图2可以看出，Δt=0.01时无论h如何取值计算均收敛，但h=1时计算精度较低，而h增大后精度得到提高；当h=10时，本文计算结果与Ghia等[14]结果吻合较好。由以上分析可知CBOS有限元法由于对流项的显式求解使得该算法受到时间步长的束缚，而当对流项采用多步格式既可降低对整体时间步长的要求又可提高计算精度。
图2 Δt=0.01时不同对流子时间步下速度沿中线的分布
Fig.2Velocity along lines through geometric center at the different h when Δt=0.01
3 带横隔板圆柱绕流特性分析
3.1 计算区域和边界条件
图3给出了带横隔板圆柱结构尺度示意图及其计算网格剖分图，LS为横隔板长度，R为圆柱半径。计算中取LS=2R(带与圆柱直径等长的横隔板)，横隔板厚度为
。取L=1 m，U=1 m/s， ρ=1.0 kg/m3，调整动力黏度得到不同的雷诺数。入口处指定沿水平方向的均匀来流u=1，垂向速度v=0；侧壁、圆柱表面和横隔板表面为固壁边界条件；出口为自由出流边界，在该位置相对压力为0。对流子时间步取h=10。
图3 带横隔板圆柱绕流结构及其计算网格剖分示意
Fig.3Size diagram and sketch of computational grid for a circular cylinder with splitter plate
3.2 结果分析
3.2.1 流场参数特性分析
表1列出了本文计算的圆柱绕流(LS =0)和带横隔板圆柱绕流平均阻力系数
、升力系数CL和斯特劳哈数St。从表1可知，在圆柱背流侧增加横隔板，当Re=100时，平均阻力系数
由1.394减小到1.223，减少了约12.3%；升力系数由0.332减少为0.076，减小了约77.1%；St由0.167减小为0.138，减少了约17.37%；在Hwang等[19]的研究中分别减少了约12.7%、76.8%和17.96%。当Re=200时，平均阻力
由1.339减小到1.0624，减少了约20.6%；圆柱所受升力幅值由0.614减少为0.163，减小了约73.5%；St由0.194减小为0.155，减少了约20.1%。计算结果表明二维层流状态下在圆柱尾迹区内放置与圆柱直径等长的横隔板使得平均阻力系数、升力系数幅值和斯特劳哈数都会减小，但都不会减为0。而段志强[4]采用二维层流模型研究Re=200的带横隔板圆柱绕流，发现横隔板长度与圆柱直径相等时升力系数幅值会减小为0。
表1 Re=100和Re=200时本文计算结果
Table 1Flow parameters for flow field at Re=100 and Re=200
参数不带横隔板(LS=0)带横隔板Re=100Re=200Re=100Re=200本文结果Han等[15]Xu和Wang[16]本文结果Ding等[17]Harichandan和Roy[18]本文结果本文结果CD1.3941.3791.4231.3391.3481.3301.2321.062CL±0.332±0.334±0.340±0.614±0.659±0.600±0.076±0.163St0.1670.1680.1710.1940.1960.1920.1380.155
3.2.2 绕流场结果特性分析
图4和图5分别给出了Re=200时圆柱绕流和带横隔板圆柱绕流半个周期内的涡量图。图4和图5对比可见：由于横隔板的存在，圆柱壁面附近的自由剪切层不能相互作用，剪切层只能沿着流动方向向下游输送，直至到达横隔板末端后，上下剪切层才开始相互作用。圆柱尾迹区的旋涡变长，相应涡脱周期变长，St变小，说明在圆柱尾部添加横隔板能有效地抑制圆柱背流面旋涡脱落。
图4 Re=200时圆柱绕流半个周期内涡量
Fig.4Vortices distribution for a circular cylinder at Re=200 during half cycle
图5 Re=200时带横隔板圆柱绕流半个周期内涡量
Fig.5Vortices distribution for a circular cylinder with splitter plate at Re=200 during half cycle
图6给出了Re=200时的圆柱绕流在大致一个涡旋脱落周期内的压力和流线的演化过程。图7给出了Re=200时带横隔板圆柱绕流一个周期内压力和流线演化。图6和图7对比可知，在两种情况下圆柱迎流面压力基本不变，而带横隔板圆柱背流面压力较不带横隔板圆柱背流面压力有着显著的增加，导致带横隔板圆柱所受平均阻力系数降低。
图6 Re=200时一个周期内压力和流线演化
Fig.6Development of pressure and streamline during a cycle at Re=200
图7 Re=200时带横隔板圆柱绕流一个周期内压力和流线演化
Fig.7Development of pressure and streamline for a circular cylinder with splitter plate during a cycle at Re=200
由图4和图5对比可见，带横隔板圆柱绕流流场尾迹区内涡量分布的对称性较不带横隔板圆柱绕流情况有了比较明显的改善，流场对称性的改善可改变圆柱上的升力系数幅值。为揭示升力幅值减小且不减为零的原因，图8给出了Re=200正向升力系数达到最大时圆柱表面压力系数Cp分布。由图8可以看出，当正向升力系数达到最大时，带横隔板圆柱背流面上部(90°
(实线：不带横隔板圆柱绕流计算结果；虚线：带横隔板圆柱绕流计算结果)
图8 Re=200正向升力系数达到最大时圆柱表面压力系数分布
Fig.8Distribution of surface pressure coefficient when the maximum of the positive lift coefficient is reached at Re=200
4 结 论
(1)CBOS有限元法中显式求解对流项在时域离散化计算中是条件稳定的，把多步格式引入到特征线算子分裂有限元法中，建立了基于多步格式的特征线算子分裂有限元法，该算法中对流项时间离散采用多步格式，在每一子时间步内沿特征线展开并显式求解。方腔流数值模拟结果证实该算法既可降低对整体时间步长的要求又可提高计算精度。
(2)带与圆柱直径等长的横隔板圆柱绕流模拟结果表明：当Re=100时，平均阻力系数
减少了12.3%；升力系数幅值减小了77.1%；St减少了17.37%；当Re=200时，3个流场参数分别减少20.6%、73.5%和20.1%。
(3)通过对比Re=200时有无横隔板绕流流场和圆柱表面压力变化表明：圆柱尾迹区放置横隔板后，剪切层只能沿着流动方向向下游输送，直至到达横隔板末端后，上下剪切层才开始相互作用，旋涡变长，旋涡脱落受到抑制，涡脱周期变长，St变小；带横隔板圆柱背流面压力较不带横隔板圆柱背流面压力有着显著的增加而迎流面压力无明显变化，导致带横隔板圆柱所受平均阻力系数降低；横隔板的加入使得圆柱上下表面压力系数差大幅度减小但仍存在差距，所以带横隔板圆柱所受升力系数幅值也大幅度降低但不会减少至0。
参考文献：
[1]吴剑, 齐鄂荣, 李炜, 等. 应用PIV系统研究横流中近壁水平圆柱绕流旋涡特性[J]. 水科学进展, 2005, 16(5): 628-633. (WU J, QI E R, LI W, et al. Research on vortexes of cross flow past a horizontal circular cylinder by PIV system [J]. Advances in Water Science, 2005, 16(5): 628-633. (in Chinese))
[2]GERRARD J H. The mechanics of the formation region of vortices behind bluff bodies [J]. Journal of Fluid Mechanics, 1966, 25(2): 401-413.
[3]KWON K, CHOI H. Control of laminar vortex shedding behind a circular cylinder using splitter plates[J]. Physics of Fluids, 1996, 8(2): 479-486.
[4]段志强. 低雷诺数下尾部隔板影响圆柱绕流的数值研究[D]. 重庆: 重庆大学, 2012: 16-25. (DUAN Z Q. An numerical investigation on fluid flow past a circular cylinder with splitter plate at low Reynolds number [D]. Chongqing: Chongqing University, 2012: 16-25. (in Chinese))
[5]KAWAI H. Discrete vortex simulation for flow around a circular cylinder with a splitter plate[J]. Journal of Wind Engineering and Industrial Areodynamics, 1990, 33(1/2): 153-160.
[6]SATAPRAHM C, LUADSONG A. The meshless local Petrov-Galerkin method for simulating unsteady incompressible fluid flow [J]. Journal of the Egyptian Mathematical Society, 2014, 22(3): 501-510.
[7]KUMAR B V R, KUMAR S. Convergence of three-step Taylor Galerkin finite element scheme based monotone Schwarz iterative method for singularly perturbed differential-difference equation [J]. Numerical Functional Analysis and Optimization, 2015, 36(8): 1029-1045.
[8]BARBOSA A G, DUTRA E G, DUDA F P. A sensitivity analysis on the parameter on the GLS for a second-gradient theory of incompressible flow [J]. International Journal for Numerical Methods in Fluids, 2015, 78(10): 638-655.
[9]TAN L, ZHU B S, WANG Y C, et al. Turbulent flow simulation using large eddy simulation combined with characteristic-based split scheme [J]. Computers& Fluids, 2014, 94(1): 161-172.
[10]ZIENKIEWICZ O C, CODINA R. A general algorithm for compressible and incompressible-flow: I: the split characteristic based scheme [J]. International Journal for Numerical Methods in Fluids, 1995, 20(8/9): 869-885.
[11]WANG D G, WANG H J, XIONG J H, et al. Characteristic-based operator-splitting finite element method for Navier-Stokes equations [J]. Science China (Technological Sciences), 2011, 54(8): 2157-2166.
[12]GLOWINSKI R, PIRONNEAU O. Finite element method for the Navier-Stokes equations [J]. Annual Review of Fluid Mechanics, 1992, 24(1): 167-204.
[13]韩先洪, 李锡夔. 高黏性流动有限元模拟的迭代稳定分步算法 [J]. 力学学报, 2006, 38(1): 16-24. (HAN X H, LI X K. An iterative stabilized fractional step algorithm for finite element analysis in high-viscosity fluid flows[J]. Chinese Journal of Theoretical and Applied Mechanics, 2006, 38(1): 16-24. (in Chinese))
[14]GHIA U, GHIA K N, SHIN C T. High-Re solutions for incompressible flow using the Navier-Stokes equations and a multigrid method [J]. Journal of Computational Physics, 1982, 48(3): 387-411.
[15]HAN Z L, ZHOU D, TU J H. Laminar flow patterns around three side-by-side arranged circular cylinders using semi-implicit three-step Taylor-characteristic-based-split(3-TCBS) algorithm [J]. Engineering Applications of Computational Fluid Mechanics, 2013, 7(1): 1-12.
[16]XU S, WANG Z J. An immersed interface method for simulating the interaction of a fluid with moving boundaries [J]. Journal of Computational Physics, 2006, 216(2): 454-493.
[17]DING H, SHU C, YEO K S, et al. Numerical simulation of flows around two circular cylinders by mesh-free least square-based finite difference methods [J]. International Journal for Numerical Methods in Fluids, 2007, 53(2): 305-332.
[18]HARICHANDAN A B, ROY A. Numerical investigation of low Reynolds number flow past two and three circular cylinders using unstructured grid CFR scheme [J]. International Journal of Heat and Fluid Flow, 2010, 31(2): 154-171.
[19]HWANG J Y, YANG K S, SUN S H. Reduction of flow-induced forces on a circular cylinder using a detached splitter plate [J]. Physics of Fluids, 2003, 15(8): 2433-2436.
*The study is financially supported by the National Natural Science Foundation of China (No.51349011).
Numerical simulation of the characteristics of flow pastcircular cylinder with splitter plate
SHUI Qingxiang, WANG Daguo, WANG Wei
(School of Environment and Resource, Southwest University of Science and Technology, Mianyang 621010, China)
Abstract:A multistep characteristic-based operator-splitting finite element method is proposed to determine the effects of adding a splitter plate at the back of a cylinder on the characteristics of a flow past a circular cylinder. In each time step, the N-S equations are split into diffusive and convective parts. For the convective part, temporal discretization is performed by the characteristic Galerkin method with a multistep technique. Results of the present model are in good agreement with previously published data. Thus, the present model is competitive in solving incompressible viscous flows in terms of accuracy and efficiency. Furthermore, the flow past a circular cylinder with a splitter plate was simulated by the present model. The flow fields and pressure on the cylinder surface with and without a splitter plate were compared. Findings indicated that the plate can suppress vortex shedding, increase the pressure on the leeward of the cylinder, and reduce the pressure difference between the upper and lower surfaces.
Key words:multistep-scheme; finite element method; driven cavity flow; flow past a circular cylinder; splitter plate
DOI：10.14042/j.cnki.32.1309.2016.04.013
收稿日期：2015-09-28;网络出版时间：2016-06-22
基金项目：国家自然科学基金资助项目(51349011)；四川省教育厅科研项目(16CZ0013)
作者简介：水庆象(1988—)，男，安徽蚌埠人，讲师，硕士，主要从事计算水力学方面研究。
E-mail:[email protected]
通信作者：王大国， E-mail:[email protected]
中图分类号：TV135;O35
文献标志码：A
文章编号：1001-6791(2016)04-0586-08
网络出版地址：http://www.cnki.net/kcms/detail/32.1309.P.20160622.2228.028.html