基于OpenFOAM编程数值模拟双圆柱绕流流动

网格
图 1 圆柱绕流计算域设置及边界条件
表 1 单圆柱绕流计算网格单元数数目与节点数目
面单元数
圆柱边界
上的
节点
数
节点 总 数
Mesh1
6550
80
6755
Mesh2
7680
96
7904
Mesh3
10600
120
10860
图 2 是Re=100 的升力和阻力系数的时间历程曲线，可以看到阻力和升力系数呈周期性振荡。表 2 给出了Re=100 计算得到的最大升力系数 C d max 、最大阻力系数 Cl max 和Strouhal数（ Str = Df / U ）。 计算结果与Schafer et al (1996) [4]给出的标准解比较也吻合的很好。图 3 为Re=100 时的涡量图。
2.2 串列双圆柱绕流
数值模拟二维串列双圆柱绕流场，重点研究了低雷诺数（ Re ≤ 500 ）下圆柱间距为 1.5D，2D，
3D，4D，5D 等几种情况下圆柱间距对流场的影响。 选取矩形计算域划分计算网格，对于串列双圆柱，入口边界距上游圆柱中心 10D，出口边界距
下游圆柱中心 20D，两侧边界距水平中线距离为 10D，上游圆柱中心位于坐标（0，0）；两圆柱中心 间距为L。入口边界条件和两侧边界条件设置为水平均匀来流U ，垂直流向速度V = 0 ；出口边界条 件设置为零法向梯度： gradϕ = 0 , (ϕ = p，U ) 。
== -fvc::grad(p)
);
代码中每一行分别与式（4）中的每一项相对应；rho → ρ ；mu → μ ；phi → φ = ρU 。
利用OpenFOAM的程序库提供的求解器和数值离散格式求解N-S方程，采用有限体积法和一次预
测两次修正的PISO算法进行计算。对N-S方程中的对流项
∇
⋅
vv (ρUU )
4
过程中，时间步长由CFL条件确定，取 Δt ＝0.001s。
(a) 全局网格
(b) 局部网格
图 4 串列双圆柱计算域网格
计算中，通过记录每个时间步的受力系数，可以得到圆柱受力的时间历程曲线。当L/D=1.5 时受
力曲线很快的达到稳定状态，而当L/D=2.0 时，受力曲线经过十几秒的不稳定振荡之后才趋于稳定状
1 数学模型与数值方法
对于不可压缩粘性牛顿流体的控制方程为Navier-Stokes方程。在直角坐标系下，连续性方程和
动量方程分别为：
∂ui = 0 ∂xi
（1）
ρ ∂ui + ρ ∂(u jui ) − ρ ∂ (ν ∂ui ) = − ∂p
∂t
∂x j
∂x j ∂x j
∂xi
（2）
式中: i, j =1，2，3； ui ,u j 为速度分量； p 为压力； ρ 为流体的密度；ν 为流体的运动粘性系数。
基于 OpenFOAM 编程数值模拟双圆柱绕流流动
曹洪建 查晶晶 万德成
（上海交通大学 船舶海洋与建筑工程学院 海洋工程国家重点实验室，上海 200240）
摘 要：本文利用开源代码 OpenFOAM 编制的 CFD 程序数值模拟双圆柱绕流流动问题，分析讨论两圆柱在串行排 列和并行排列情况下，间距的改变对绕流流场和各圆柱受力情况的影响。数值模拟结果表，基于 OpenFOAM 的开源 程序库编制的 CFD 程序，可以很好和灵活处理复杂的流动问题，为开发 CFD 新格式和新方法提供了基础。 关键词：OpenFOAM，双圆柱，串列，并列，CFD
选取网格数量不同的三种四边形规则结构网格进行计算。三种网格的数量和圆柱边界上的节点 数在表 1 中给出。
2
本文数值模拟Re=20，100 时的单圆柱绕流场。Re=20 的流场稳定时，圆柱的后面附着一对固定 的漩涡，不发生漩涡脱落。当Re =100 时，圆柱后的尾涡交替脱落，同时又生成新的尾涡，形成涡街。
网格
图 2 单圆柱绕流的阻力系数、升力系数的时间历程曲线（Re=100）
表 2 Re=100 时数值计算结果及与标准解的相对误差
最大阻 力 系
相对 误 差
最大升 力 系
相对 误 差
数
数
Error
E
Cd max
1
Cl max
r
3
Strouha l 数
Str
Mes h 1
Mes h 2
Mes h 3
标准 解
态。 L/D=3.0 时，受力曲线也很快达到稳定，而当L/D ≥ 4.0 时，受力曲线达到稳定状态的经历的时
间延长，在 L/D=5.0 时，受力曲线在时间超过 25 秒之后才逐渐趋于稳定状态。图 5 为Re=200，L/D=2.0
和L/D=4.0 时计算得到的阻力和升力系数的时间历程曲线。
(a) L / D=2.0 5
网格的质量和数量将影响数值模拟的精度和计算速度。对于双圆柱计算域，本文采用混合结构 网格分区进行划分。以串列圆柱间距L=3D的网格为例，整个计算域的网格划分如图 4(a)所示；在圆 柱周围采用贴体的四边形结构网格，如图 4(b)所示，圆柱壁面上均匀布置 120 个节点，靠近圆柱壁 面最近一层的网格节点到圆柱壁面距离为 0.021D，两圆柱之间以及圆柱的下游采用三角形非结构网 格；在远离圆柱的上游和两侧部分采用四边形结构网格，通过减少网格的数量减少计算时间。计算
采用四阶精度的高斯积分离散格
式（Gauss
cubic）；拉普拉斯项
∇
⋅
v (μ∇U )
采用高斯积分的线性修正离散格式（ted）；
压力梯度项采用高斯积分的线性离散格式（Gauss linear）；时间上采用欧拉隐式离散格式。
2 数值计算与结果讨论
本文研究的不可压缩粘性流体的流体密度 ρ = 1.0kg/m3 ，流体运动粘性系数ν = 10-3 m2/s，所有 算例中的圆柱直径 D = 0.1m ,U 表示平均速度，定义雷诺数 Re =UD /ν 。
表 3 串列双圆柱绕流的数值计算结果
L/D
雷诺数 Re
平均阻力系数 Cd1
升力系数 Cl1
平均阻力系数 Cd2
阻力系数 Cl2
Strouhal 数 St
1.5
40
100
200
500
1.582 1.240 1.116 1.057
-0.00013
± 0.018 ± 0.027 ± 0.14
0.122 -0.091 -0.207 -0.307
1.19e-4
± 0.053 ± 0.068 ± 0.398
引言
圆柱绕流是流体力学的经典问题之一，广泛地存在于海洋工程领域中，随着对海洋海底资源的 开发，人们越来越关注海洋平台立管和海底管线等在水流作用下的受力和振动所带来的安全问题， 在海洋工程中，许多结构物是由多个圆柱管件构成，因此对多圆柱绕流的研究也非常重要，具有更 重要的工程应用意义。
国内外众多学者作了大量的研究工作研究均匀来流中固定单圆柱的二维绕流，对单圆柱绕流问 题已有较多的结果。对双圆柱绕流的研究可以更深刻的了解圆柱间的干扰和由圆柱泻涡引起的振动， 因此，国内外学者对多圆柱绕流问题也进行了研究。Zdravkovich[1] 曾对均匀流中并列、串列和斜向 交错布置等情况下的双圆柱绕流问题进行了研究，研究表明在两个圆柱串列情况下，存在三种不同 流态：当两圆柱中心距离 L<(1.2D-1.8D)时，串列圆柱表现为单一钝体绕流，从上游圆柱脱离的剪切 层不会在下游圆柱上附着；当(1.2-1.8)<L/D<(3.4-3.8)时，从上游圆柱脱离的剪切层会附着在下游圆 柱 上 ， 尾 涡 的 脱 落 仅 有 下 游 圆 柱 产 生 ； 当 L>4D 时 ， 两 个 圆 柱 都 会 发 生 尾 涡 脱 落 现 象 。 Bearman&Wadcock[2]和Williamson[3]对并列双圆柱绕流进行了研究，结果表明：当L/D<2.2时，两圆 柱只形成一个涡街，并发生共振情况。当L/D>2.2时,形成两个相位相反的涡街，当L/D>5，随着L的 增大，绕流情况逐渐与单个圆柱的绕流情况相似。
3.2057 2
3.2086 4
3.2127 9
3.2200 ～ 3. 24 00
0.752 %
0.661 %
0.533 %
――
0.9561 71
0.9690 63
0.9834 06
0.9900 ～ 1. 01 00
r o r 2 4.329 %
3.037 %
1.659 %
――
0.2971
0.2966
首先，对二维单圆柱绕流场进行数值模拟，并与标准解进行比较，验证基于OpenFOAM进行圆 柱绕流场数值模拟的有效性和准确度。
2.1 单圆柱绕流 数值模拟单圆柱绕流场是为了检验应用OpenFOAM程序求解Navier-Stokes方程的有效性和计算
精度。本文采用Schafer et al (1996) [4]提供的计算模型，二维计算域的设置和边界条件如图 1 所示， 计算域长度 L = 2.2m ,宽度 H = 0.41m ，圆柱的直径 D = 0.1m ，圆柱中心坐标为 (0 .2,0 .2 ) 。入口为抛 物 线 型 来 流 U (0, y) = 4U m y( H − y) / H 2 ， V = 0 , 其 中 U m = U (0, H / 2) , 平 均 速 度 U = 2U (0, H / 2) / 3 = 2U m / 3 ；圆柱表面和两侧边界为无滑移壁面条件U = V = 0 。
0.2970 0.2950
～ 0. 30 50
应用OpenFOAM进行数值模拟二维单圆柱绕流，计算得到的圆柱的阻力、升力和Strouhal数等结 果与Schafer et al (1996) [4]的标准解相吻合，验证了本文应用OpenFOAM 求解器数值模拟圆柱绕流场 的有效性和准确性。