自动网格体系在柱体绕流大涡模拟中的适用性评估

摘要

为了评估基于snappyHexMesh方法生成的自动网格体系在二维柱体绕流大涡模拟中的适用性，比较了该自动网格体系与人工网格体系对于Re为3 900圆柱绕流和Re为22 000方柱绕流的数值模拟结果。通过设置合理的计算域以及数值格式，采用snappyHexMesh自动网格以及人工网格的算例都表现出良好的数值稳定性。将不同网格体系的数值模拟结果与物理试验结果进行对比，结果表明，采用snappyHexMesh网格可以提高数值求解效率；圆柱绕流对网格体系的变化比较敏感，不同密度的snappyHexMesh网格会显著影响圆柱气动力特征以及尾流区域的流场结果；snappyHexMesh网格体系可以准确预测方柱绕流，在方柱绕流大涡模拟中具有相较于圆柱绕流更好的适用性。

关键词

计算流体力学; 自动网格体系; 柱体绕流; 大涡模拟; 气动力特征; 流场结构

采用计算流体力学（computational fluid dynamics， CFD）的方法解决实际工程问题时，CFD数值计算结果的准确性会受到诸如网格划分体系、数值计算格式、求解精度等因素的影响而发生改变^［1］。其中，网格划分体系作为数值计算的重要部分之一，其质量的高低会直接决定算例的计算效率、计算稳定性及计算结果的可信程度。风工程中常用的大涡模拟（large eddy simulation， LES）方法是一种对网格体系质量较为敏感的CFD方法^［2］，针对网格体系在大涡模拟中的适用性，研究者开展了大量数值模拟研究，Wang等^［3］比较了六面体、四面体和多面体网格体系在建筑物周围区域边界层风场数值模拟中的适用性；Mukha等^［4］比较了底面为三角形，四边形和多边形的棱柱网格对壁面区域LES计算精度的影响。由于网格质量没有唯一的评价标准，当采用一种新的网格体系时，对该网格体系在LES计算中的适用性开展研究具有十分重要的工程意义。

按照不同的划分方法可以把网格体系分为两类：人工网格体系和自动网格体系。人工网格体系是指在计算开始之前，根据模拟对象和计算目的人为地规划网格分布，由于其质量可控且适用性较强，被广泛用于CFD研究，例如：周强等^［5］采用人工划分方式得到了过渡均匀的六面体网格体系来研究中高雷诺数下的圆柱绕流特性；王蒙恩等^［6］利用人工方法对真实龙卷风模拟装置建立数值模型并开展大涡模拟。然而在较为复杂的结构绕流计算中，采用人工划分网格的方法技术难度较高。为了提高实际工程计算中的网格划分效率，可以采用基于一定规则设计的自动网格划分方法对计算域进行离散，获得满足计算需要的自动网格体系。国内外学者都基于自动网格体系开展大量研究，例如；Ricci等^［7］在模拟边界层风场与低矮建筑物的相互作用时采用自动网格划分方法对计算域进行离散；Du等^［8］提出了一种适用于复杂城市区域的高效的网格自动生成方法，将其应用于香港工业大学的行人高度处风环境数值模拟。目前已经有许多开源或商用的自动网格划分脚本或软件，比较著名的有商业软件Fluent、Pointwise等。流体计算开源平台OpenFOAM提供了一种三维网格自动划分程序snappyHexMesh（以下简称Hex），该程序能够基于输入的参数自动生成分布均匀的贴体六面体网格体系，其网格生成速度快且支持并行计算，非常适用于大尺度密网格的复杂区域流场计算，目前已经有许多成功的应用，例如：Cavar等^［9］基于Hex程序对复杂山区地形进行建模并对该地区的平均和脉动风场进行预测；Rakai等^［10］采用不同网格划分体系对某一理想化城市街区的风环境及污染物扩散开展数值模拟，结表明Hex网格可以取得较好数值结果和计算稳定性。

为了进一步研究在基于LES方法的数值模拟中Hex网格对流场结果的影响，本文从基础的流体绕流现象着手，分别对雷诺数Re为3 900的均匀来流下圆柱绕流和Re为22 000的均匀来流下方柱绕流进行不同网格体系下的LES模拟，通过比较采用人工网格体系和采用Hex网格体系得到的柱体气动力特征、计算域内流场的瞬时结果和平均结果，定量分析Hex网格在平面和曲面绕流数值模拟中的适用性，为Hex网格体系在一般工程问题中的应用提供参考。本文首先介绍了研究目的和研究内容；之后在第一部分介绍了数值模拟的前期准备工作；第二部分展示了圆柱和方柱绕流计算的不同网格体系，包括人工划分的平滑过渡网格和自动划分的Hex网格；第三部分详细比较了各组数值模拟结果的计算效率、柱体气动力特征及流场特征；第四部分对全文进行了总结。

1 数值模拟方法及计算域设置

1.1　大涡模拟方法的控制方程

大涡模拟是由Smagorinsky^［11］提出的一种湍流数值模拟方法，其基本思想是通过指定空间滤波函数将湍流分解为可解尺度湍流和不可解尺度湍流脉动，对可解尺度湍流脉动进行直接数值模拟，对不可解尺度湍流脉动采用亚格子模型来描述。经过滤波操作后，大涡模拟（LES）方法的连续性和N-S方程为

$$\frac{\partial \rho }{\partial t}+\frac{\partial \left(\rho {\overline{u}}_i\right)}{\partial x_i}=\mathrm{0}$$

（1）

$$\frac{\partial \left(\rho {\overline{u}}_i\right)}{\partial t}+\frac{\partial \left(\rho {\overline{u}}_i{\overline{u}}_j\right)}{\partial x_j}=\frac{\partial \overline{p}}{\partial x_i}+\frac{\partial }{\partial x_j}\left(\mu \frac{\partial {\overline{u}}_i}{\partial x_j}\right)-\frac{\partial {\tau }_{ij}}{\partial x_j}$$

（2）

式（1）—（2）中：$\rho$为流体密度；$u_i$、$u_{\mathrm{j}}$为速度分量；$p$为压力；$\mu$为流体运动粘度；其中带有上划线的量为滤波后的场变量；${\tau }_{ij}=\overline{u_i{u}_j}-{\overline{u}}_i{\overline{u}}_j$被定义为亚格子应力。常用的Smagorinsky-Lilly^［12］亚格子模型因其形式简单，数值运算稳定，成为一种广泛应用的亚格子模型。在Smagorinsky-Lilly亚格子模型中：

$${\tau }_{ij}-\frac{\mathrm{1}}{\mathrm{3}}{\delta }_{ij}{\tau }_{kk}=-\mathrm{2}{\left(C_{\mathrm{s}}\overline{\Delta }\right)}^{\mathrm{2}}\left|{\overline{S}}_{ij}\right|{\overline{S}}_{ij}$$

（3）

式中：Smagorinsky常数$C_{\mathrm{s}}=\mathrm{0.1}$；$\overline{\Delta }$为网格滤波尺寸，本文采用立方体滤波方法：$\overline{\Delta }=\sqrt[\mathrm{3}]{∆x∆y∆z}$；${\overline{S}}_{ij}=\frac{\mathrm{1}}{\mathrm{2}}\left(\frac{\partial {\overline{u}}_i}{\partial x_j}+\frac{\partial {\overline{u}}_j}{\partial x_i}\right)$，$\left|\overline{S}\right|=\sqrt[]{\mathrm{2}{\overline{S}}_{ij}{\overline{S}}_{ji}}$。为了提高近壁面处的湍流模拟结果，本文采用Van Direst壁面函数对标准Smagorinsky-Lilly模型进行修正。

1.2　计算域设置及求解格式

本文共模拟两类柱体绕流现象，分别为均匀来流下的亚临界雷诺数圆柱绕流以及高雷诺数方柱绕流，每类柱体绕流包含三组不同算例，各组算例除了网格体系存在区别外，其计算域尺寸（除了展向高度）、边界条件（除了入流速度）、求解方法和计算格式完全一致，分别介绍如下。

对于计算域尺寸，本文各组算例采用立方体计算域，顺流向长度30 D（其中D为圆柱断面直径或方柱断面宽度），横流向长度17 D，圆柱或方柱模型的上游来流区域长度为10 D，下游尾流区域长度为20 D，圆柱或方柱的几何中心距离计算域的上下边界各8.5 D。对于圆柱模型其计算域展向高度为$\mathrm{\pi }$D，展向均匀划分48层网格；对于方柱其计算域展向高度为6 D，展向均匀划分60层网格。计算域示意如图1所示，数值模型的边界条件如表1所示。

图1  计算域示意图及边界条件

Fig. 1  Computational domain and boundary conditions

表1  模型边界条件汇总

Tab. 1  Summary of boundary conditions

边界名	速度边界条件	压力边界条件
入流边界	定值	零梯度
出流边界	零梯度	定值
展向边界	周期边界条件
前后边界	对称边界条件
柱体表面	无滑移边界条件

在求解格式方面，本文采用了二阶精度的数值算法，对流体方程的时间项差分采用backward二阶隐式格式，压力梯度项采用基于高斯积分的二阶线性差分格式，对流项采用线性迎风二阶差分格式，扩散项采用基于高斯积分的二阶线性修正差分格式。对压力速度耦合方程采用隐式压力分裂算子（PISO）算法求解，并在每个求解时间步内进行2次非正交修正来保证计算稳定性。

2 网格划分体系

2.1　圆柱绕流网格体系

为了比较自动网格划分体系对亚临界雷诺数（Re为3 900）圆柱绕流计算结果的影响，本文设置了三个算例：第一个算例采用人工方法划分网格，得到了过渡平滑的结构化网格，网格数量约为300万；第二个和第三个算例采用程序自动划分的方式得到了Hex网格体系，其中第二个算例的Hex网格整体密度与人工网格体系接近，网格数量约为550万，而第三个算例采用网格相对稀疏但网格数量与人工网格体系接近的Hex网格体系，约为350万。为了保证LES方法在圆柱表面边界层内计算的准确性，三组网格体系下的圆柱表面网格高度均为0.002 D，使得计算过程中的$y^{+}$最大值小于1。三组网格体系的整体及局部示意如图2所示，网格尺寸沿流向的过渡情况如图3所示。

图2  圆柱绕流计算的不同网格体系示意图

Fig. 2  Different grid systems for simulation of flow around circular cylinder

图3  圆柱绕流计算不同网格体系的尺寸沿流向过渡特征

Fig. 3  Change of grid size along flow direction for simulation of flow around circular cylinder

2.2　方柱绕流网格体系

对于方柱绕流的计算，本文同样地设置三组算例来比较自动网格划分体系对高雷诺数（Re为22 000）方柱绕流计算结果的影响，其网格数量分别约为330万，800万和410万。三组网格体系下的方柱表面网格高度分别为0.002 D、0.0015 D和0.003 D，使得计算过程中的$y^{+}$最大值小于5。用于方柱绕流计算的三组网格体系的整体及局部示意如图4所示，网格尺寸沿流向的过渡情况如图5所示。

图4  方柱绕流计算的不同网格体系示意图

Fig. 4  Different grid systems for simulation of flow around square cylinder

图5  方柱绕流计算不同网格体系的尺寸沿流向过渡特征

Fig. 5  Change of grid size along flow direction for simulation of flow around square cylinder

3 计算结果分析

3.1　数值求解效率

本文的各组绕流算例都采基于 OpenFOAM V2.3.1求解，采用64核服务器并行计算，并通过调整算例的时间步长使得计算全程的库朗数$C_{\mathrm{F}\mathrm{L}}\approx$1，提取计算稳定后40 s的结果进行分析。各算例的计算总时长$T$、总时间步$n$，总网格单元数$k$以及每百万网格每时间步平均求解时长（$T/nk$，用来表征数值求解效率）如表2所示。对于两类不同的绕流计算，采用人工网格体系所需的百万网格平均每步求解时长都高于采用Hex网格的时长，即采用Hex网格体系的数值求解效率较高；对比不同密度的Hex网格的结果发现降低网格密度可以明显提高数值求解效率。结果表明，由于Hex网格体系具有优异的正交特性，当采用Hex网格体系对平直表面进行离散时可以改善由于人工网格正交性较差导致的完成计算所需的总时间步增加的现象，进而降低总计算时长；同时Hex网格优异的网格正交性也可以提高数值求解的效率；但是Hex网格对弯曲边界的适应性较差。

表2  各算例计算时长及效率比较

Tab. 2  Comparison of total solving time and efficiency among cases

算例		计算总时长T/s	总时间步n	总网格数k/百万	每百万网格平均每步求解时长/s
圆柱	人工网格	63 925	50 000	3.0	0.43
	密集Hex网格	221 927	100 000	5.5	0.40
	稀疏Hex网格	126 961	100 000	3.5	0.36
方柱	人工网格	301 753	200 000	3.3	0.46
	密集Hex网格	399 923	133 333	8.0	0.37
	稀疏Hex网格	84 751	80 000	4.1	0.26

3.2　气动力特征

对于亚临界雷诺数的圆柱绕流，图6展示了三组网格下圆柱的阻力系数$C_{\mathrm{D}}$和升力系数$C_{\mathrm{L}}$的时程结果，其中力系数定义为$C_x=\mathrm{2}{F}_x/(U_{\mathrm{\infty }}^{\mathrm{2}}D)$，$F_x$为柱体受到的力的大小，下标表示升力或阻力，$U_{\mathrm{\infty }}$为入流边界的顺流向流速。网格体系对圆柱的升力和阻力时程有明显影响，特别是采用稀疏Hex网格会显著加剧圆柱的阻力和升力脉动特性；图7展示了圆柱表面的平均压力系数${\overline{C}}_{\mathrm{p}}$的分布情况，其中压力系数定义为$C_{\mathrm{p}}=\mathrm{2}(p-p_{\mathrm{\infty }})/U_{\mathrm{\infty }}^{\mathrm{2}}$，$p_{\mathrm{\infty }}$为入流边界角点处压力。从结果来看，采用密集Hex网格获得了与试验值最接近的结果，而稀疏Hex网格略微高估了圆柱背风区的负压绝对值。表3给出了不同网格体系得到的基于时间积分的圆柱气动力特征参数，包括分离角${\theta }_{\mathrm{s}\mathrm{e}\mathrm{p}}$，斯托洛哈数St，阻力系数平均值${\overline{C}}_{\mathrm{D}}$，升力系数脉动值$C_{\mathrm{L}-\mathrm{R}\mathrm{M}\mathrm{S}}$，平均背压${\overline{C}}_{\mathrm{p}-\mathrm{b}\mathrm{a}\mathrm{c}\mathrm{k}}$，回转长度$L_{\mathrm{r}}$（定义为圆柱体后边缘到中心线上流速首次由负转为正对应位置的距离） 和柱后方最小速度与平均速度的比$U_{\mathrm{m}\mathrm{i}\mathrm{n}}/U_{\mathrm{\infty }}$，结果表明，本文的数值模拟结果与既有结果较为符合，证明了本文计算设置的合理性以及Hex网格的适用性。

图6  圆柱绕流阻力系数和升力系数时程结果比较

Fig. 6  Comparison of time-varying drag and lift coefficients for simulation of flow around circular cylinder

图7  圆柱表面平均压力系数分布

Fig. 7  Distribution of mean pressure coefficient on circular cylinder

表3  圆柱绕流气动力特征参数比较

Tab. 3  Comparison of aerodynamic characteristic parameters for flow around circular cylinder

算例	${\theta }_{\mathrm{s}\mathrm{e}\mathrm{p}}$	St	${\overline{C}}_{\mathrm{D}}$	$C_{\mathrm{L}-\mathrm{R}\mathrm{M}\mathrm{S}}$	$-{\overline{C}}_{p-\mathrm{b}\mathrm{a}\mathrm{c}\mathrm{k}}$	$L_{\mathrm{r}}$	$-U_{\mathrm{m}\mathrm{i}\mathrm{n}}/U_{\mathrm{\infty }}$
人工网格	87.75	0.212	0.987	0.107	0.943	1.49	0.22
密集Hex网格	86.80	0.213	0.970	0.081	0.877	1.71	0.31
稀疏Hex网格	88.21	0.210	1.062	0.240	0.987	1.31	0.30
Ong和Wallace［14］（实验）	85±2	0.215±0.005	0.99±0.05	—	0.9±0.05	1.33±0.05	0.24±0.01
Lourenco和Shih［15］（实验）	86˚	0.22	0.99	—	—	1.19	0.24
Mani等 ［16］（LES）	86˚	0.21	0.99	—	0.86	—	0.33
Wornom等［17］（LES）	89˚	0.21	0.99	0.11	0.88	1.45	—

与圆柱类似的，对于方柱绕流，图8展示了三组网格下方柱的阻力系数$C_{\mathrm{D}}$和升力系数$C_{\mathrm{L}}$的时程结果，不同网格体系对方柱计算的影响较小；图9和图10分别展示了方柱表面的平均压力系数${\overline{C}}_{\mathrm{p}}$及脉动压力系数${\overline{C^{\text{'}}}}_{\mathrm{p}}$的分布情况，由于方柱绕流的漩涡分离点固定，其表面压力分布被明显地分成4部分，对于迎风面，三组不同网格体系的算例得到的表面压力结果一致且与Nishimura和Taniike^［18］以及Oka等^［19］的结果吻合；对于两个侧风面，三组计算的平均风压绝对值偏高，而脉动风压结果与试验结果比较符合，三组网格对侧风面表面风压的影响不明显；而对于背风区，三组网格结果的表面平均风压绝对值和风压脉动绝对值都偏低，采用稀疏Hex网格使得风压结果更接近试验值。表4汇总了不同网格体系得到的基于时间积分的方柱气动力特征参数，各项气动力统计值与既有结果比较符合，但是其回转长度结果都偏低，产生偏差的原因可能是由于不同文献中对回转长度的定义不一致。

图8  方柱绕流阻力系数和升力系数时程结果比较

Fig. 8  Comparation of time-varying drag and lift coefficients for simulation of flow around square cylinder

图9  方柱表面平均压力系数分布

Fig. 9  Distribution of mean pressure coefficient on square cylinder

图10  方柱表面脉动压力系数分布

Fig. 10  Distribution of pressure fluctuation coefficient on square cylinder

表4  方柱绕流气动力特征参数比较

Tab. 4  Comparison of aerodynamic characteristic parameters for flow around square cylinder

算例	Re	St	${\overline{C}}_{\mathrm{D}}$	$C_{\mathrm{D}-\mathrm{R}\mathrm{M}\mathrm{S}}$	$C_{\mathrm{L}-\mathrm{R}\mathrm{M}\mathrm{S}}$	$-{\overline{C}}_{\mathrm{p}-\mathrm{b}\mathrm{a}\mathrm{c}\mathrm{k}}$	$L_{\mathrm{r}}$
人工网格	22 000	0.121	2.316	0.271	1.526	1.584	1.110
密集Hex网格	22 000	0.118	2.175	0.186	1.479	1.481	1.125
稀疏Hex网格	22 000	0.121	2.329	0.228	1.608	1.618	1.050
Nishimura和Taniike［18］（实验）	40 000	—	2.33	0.26	1.33	1.60	—
Norberg［20］（实验）	22 000	0.130	2.10	—	—	1.37	—
Wang等［21］ （ LES）	21 400	0.130	2.03	0.18	1.29	—	1.31
周强等［22］ （ LES）	22 000	0.122	2.27	0.29	1.40	1.61	1.37

3.3　流场特征

图11给出了圆柱和方柱绕流在不同网格体系下的柱体尾流区中心线的平均流向速度${\overline{U}}_x$与入流平均速度$U_{\mathrm{\infty }}$的比随距离的分布结果。对于圆柱绕流尾流区域内，采用稀疏Hex网格的结果与Lourenco和Shih^［15］的物理实验结果更为接近，而采用密集Hex网格的结果更接近于Parnaudeau等^［23］的试验结果，采用人工网格的结果介于二者之间；对于方柱绕流，三组网格的尾流区结果几乎一致，且与Wang等^［21］的大涡模拟结果非常接近，然而与Lyn等^［24］的试验结果有一定差距。结果表明，相比于方柱绕流，圆柱绕流尾流区域对网格体系的改变更加敏感，改变网格的划分方法或者疏密程度都会显著改变圆柱绕流的流场结果。为了更直观地比较不同网格体系对于圆柱绕流周围流场的影响，图12分别展示了三组圆柱算例的时均流线结果，人工网格成功再现了大尺度的回流结构，同时还捕捉到了附着在圆柱表面的小尺度涡旋结构^［25］；采用密集Hex网格得到的回流区域较大且回流漩涡距离圆柱较远，同时小尺度涡旋消失；而采用稀疏Hex网格得到的回流区较小且距离圆柱较近，同时圆柱表面的再附着涡旋强度减弱且范围增大。

图11  尾流中心线上平均顺流向速度分布

Fig. 11  Distribution of mean streamwise velocity in wake centerline

图12  圆柱绕流尾流区域平均流线分布

Fig. 12  Time-averaged streamlines in the wake of cylinder

为了比较柱体后方区域的横流向流场速度分布情况，图13给出了圆柱和方柱绕流在不同网格体系下的柱体尾流区不同位置处平均顺流向速度的分布结果。平均流向速度呈现“U”形分布，且并随着流动方向逐步趋于平缓，采用人工网格或者不同密度的Hex网格的顺流速度偏差仅集中在尾流区域，对圆柱和方柱周围其他区域的顺流向流速场影响较小。

图13  尾流区域不同位置的平均顺流向速度分布

Fig. 13  Mean streamwise velocity at different locations in the wake of cylinder

在绕流过程中，由于柱体两侧存在漩涡交替脱落的现象，会导致在柱体后方的流场中出现明显的顺流向脉动速度。图14分别给出了圆柱和方柱绕流在不同网格体系下的柱体尾流区不同位置处顺流向脉动速度均方值的分布结果，图中$u^{\text{'}}$为顺流向脉动速度。对于圆柱绕流，其顺流向脉动速度随网格的改变而明显不同，特别是稀疏Hex网格明显低估了尾流区的脉动强度；对于方柱绕流，三组计算的结果与Lyn等^［24］的试验结果一致，且明显优于Srinivas等^［26］的LES结果，网格体系对方柱尾流区域顺流向脉动的影响不明显。

图14  尾流区域不同位置的顺流向脉动速度均方值分布

Fig. 14  Mean square value of streamwise velocity fluctuations at different locations in the wake of cylinder

3.4　结果讨论

对于大涡模拟来说，本文采用具有二阶精度的数值格式开展计算，证明了Hex网格对于较高精度数值离散格式具有较好的适用性。然而，当采用相对低阶的数值格式进行大涡模拟时，需要对Hex网格的适用性重新展开评估。采用一阶数值离散格式时相同网格体系的圆柱绕流计算结果如图15所示，从结果来看，使用一阶精度数值格式时采用Hex网格会得到错误的圆柱表面压力系数预测结果，这说明Hex网格对低阶精度数值格式的适用性较差。从另一方面来讲，高阶精度的数值离散格式对于网格体系的依赖性较弱，而低阶离散格式对网格的依赖性相对较强，需要配合高质量的网格体系才能保证求解的正确性。然而关于不同的数值离散方法与Hex网格的依赖性的讨论仍有待后续的进一步研究。

图15  不同数值离散精度下圆柱表面平均压力系数分布

Fig. 15  Distribution of mean pressure coefficient adopting different numerical discretization accuracies

4 结论

本文基于典型二维亚临界雷诺数圆柱绕流和高雷诺数方柱绕流现象，评估一种采用snappyHexMesh（简称Hex）方法生成的自动网格体系在二维绕流大涡模拟中的适用性。通过将三组不同的网格体系下的数值模拟结果与物理试验等结果进行对比，证明了Hex网格用于二维绕流大涡模拟计算的可行性，评估了Hex网格体系的数值求解效率，并比较了网格对于柱体绕流气动力特征和流场特征的影响，主要结论如下：

（1）采用Hex网格会提高数值求解的效率；相比于平直边界，Hex网格在弯曲边界处的正交性较差，因此需要较短的时间步长来保证计算稳定性，增加计算所需的总时间。

（2）网格体系的改变会对圆柱绕流的气动力特征产生明显的影响，采用高分辨率的Hex网格可以得到与试验值最接近的表面压力结果；网格体系对方柱绕流的影响不明显。

（3）圆柱绕流的尾流区域流场特征的大涡模拟结果对于网格体系的改变比较敏感；而不同网格对方柱绕流模拟的影响较小。采用低分辨率Hex网格对尾流区域脉动风速的预测结果较差。

作者贡献声明

张宇鑫：实验思路构思，数据的分析，实验结果可视化，初稿撰写，初稿修改。

曹曙阳：方法指导，试验数据核实，初稿审阅和修改。

操金鑫：方法指导。