网格对CFD模拟结果的影响分析

刘明亮，张思青，李胜男

（1.昆明理工大学冶能学院，云南省昆明市 650093;2.云南电网有限责任公司电力科学研究院，云南省昆明市 650000）

网格对CFD模拟结果的影响分析

刘明亮1，张思青2，李胜男3

（1.昆明理工大学冶能学院，云南省昆明市 650093;2.云南电网有限责任公司电力科学研究院，云南省昆明市 650000）

在CFD模拟的过程中，高质量的网格对计算结果的影响至关重要。本文通过一个继电器内部流场的实例，描述了网格独立性检验的详细过程。通过Fluent双精度求解器，应用k-omega模型对比添加与不添加边界层时计算误差，进而探讨边界层网格对计算精度的影响。

网格无关解；网格独立性检验；k-omega模型；边界层网格划分

0 引言

作为流体动力学研究方法之一的CFD模拟法，随着理论研究和实际应用的深入以及商用软件的日渐成熟，越来越受到工程应用者的重视，甚至产生一定的依赖。这就对计算可信度的要求也越来越大，国内已有大量的文献对不同方法给出误差的理论分析[1]-[5]。

但不管理论如何完善，实际计算中的大量误差还是避免不了的，如不确定的边界条件、近似的计算模型、简化的几何模型、低质量的网格等，以上几个因素的综合，可能会导致计算结果毫无意义。通常的解决办法是配以模型试验的数据验证。然而，并不是所有的CFD计算都有可行的实验来验证，数值模拟结果的误差本身是可以接受的，不过这需要CFD软件应用者能深刻的理解计算模型的原理，每个模型的适用条件，从而可以最大程度的减少人为误差。

实际计算中，我们可以用商用软件画出质量很好的网格，但是网格质量高，并不代表网格数量足够满足计算精度的要求。找到合适计算的网格数量需要网格独立性检验，进而得到的结果叫网格无关解。本文就是在假定模型选择合理的条件下，研究网格对计算结果的影响。

1 误差理论

CFD的模拟就是用计算机数值方法，求解控制流体运动的包含初边值问题的偏微分方程。若不考虑模型误差，数值计算中代入的误差主要包括：用有穷表示无穷而产生的截断误差和由计算机的位宽有限因估值而产生的舍入误差。

CFD计算实质上是将连续的初边值问题离散化，转化为离散的初边值问题。常见的离散方法有有限元、有限体积、有限差分等方法，理论上都是网格数越多，结果越接近真值。当网格无限多且无限小时，差分方程将变成精确的理论解。

实际上：网格不可能无限制的加密。主要的问题有：首先，网格越密，则计算量越大，计算周期也越长，通常我们的计算资源总是有限的。其次，随着网格的加密，为了计算的比敛时间步长也必须减小。实际的经验告诉我们当网格尺寸及时间步长减小时，截断误差确实减小，而舍入误差却反而增加（如图1所示）。持续减小步长尺寸并不意味着能获得更加精确的结果。事实恰好相反，在很小的时间步长处，因为舍入误差的增加，使计算精度反而下降。计算机浮点运算造成的舍入误差也会增大。因此在实际应用中，我们寻找的就是计算精度与计算开销间一个平稀点，也就是达到网格无关的阈值[6]。

所谓网格无关性检验，就是验证网格数独立于计算精度。通常是找到一个满足CFD计算的最少网格数。执行时，让网格密度逐渐增加，并选择一个判断参数，观察其随网格数量增加的变化情况，当网格的增加几乎不再影响该参数的变化时，网格数量达到计算要求。

图1 网格尺寸或时间步长与离散误差、舍入误差、总误差的函数关系Fig.1 Function relation between mesh size or time step and discrete error，rounding error，total error

2 计算实例

2.1 几何模型

本算例所选计算模型是国产沈阳四兴QJ-80型瓦斯继电器，瓦斯继电器是变压器的非电量保护重要的元器件，用三维绘图软件ProE建立的几何模型，用ANSYS中的几何处理模块Design Modeler提取几何流道生成流体计算域如图2所示，相应的部件部分经过简化，标记在该图上。

图2 瓦斯继电器三维模型（上）与流体计算域三维模型（下）Fig.2 Buchholz relay three-dimensional model（on）and computational fluid domain three-dimensional model（down）

继电器的工作原理见示意简图（如图3所示），当变压器内部发生严重故障（绝缘击穿、匝间短路、铁芯事故等）时，产生强烈的瓦斯气体，变压器油箱内压力瞬时突增，产生很大的流向油枕方向的油流冲击，冲击继电器中挡板，产生挡板动作，进而触发报警信号。

图3 瓦斯继电器工作原理示意图Fig.3 Schematic diagram of the Buchholz relay working principle

本算例的目的是为精确计算挡板上油流冲击产生的阻力力矩而准备高质量网格。下面先从没有边界层的网格起，到添加边界层的网格，用ICEMCFD软件来详细说明准备计算网格的过程。

2.2 无边界层网格的生成

进行网格独立性检验首先要确定判断参数。通常不可压计算，选用进出口压差来作为判断收敛的参数。当进口的压强不变时，压差只取决于出口值，因此本算例用出口压强面积加权平均值来作为判断收敛参数。

其次要绘制初始网格，初始的网格不需要太细致。模型的入口管径为80mm，因为我们关心挡板处的结果，因此将挡板处最大网格边长设为5mm，其余部分最大网格边长定为8mm。生成的网格记为“网格1”。

下面对网格进行加密。ICEM-CFD中，功能“Scale Factor”用来改变全局的网格尺寸（体、面、线），通过将它乘以设置的参数来得到实际网格参数。即可以通过调整Scale Factor值的大小来改变网格密度。

我们选用加倍的网格加密原则，即加密后的网格数量是加密前的两倍。为达到该目的，只需另加密后的边长变为原来的，这个值约等于0.7937。第二次加密就使Scale Factor变为0.79372，以此类推，第n次加密，Scale Factor 的值就为0.7937n。

以网格1为基础，依据以上原则，生成网格2，网格3，一直到网格5，每次画好的网格再经过光顺，保证网格质量在0.4以上。用 fluent计算每个网格，观察出口相对静压值，并将相应的数据汇总于表1中，并以图示的形式展示在图4中。

表1 出口相对静压与网格数的关系Tab.1 Relationship between export relatively static and number of grid

图4 出口相对静压随网格数的变化Fig.4 Exports relatively static pressure changes with the number of grid

从图4可以看出当网格达到1337948时，出口压强平均值已经逐渐稳定。从网格4到网格5，出口压强变化很小，也就是说，当网格数量达到1337948时，再增加网格数量已经不会显著提高计算精度，进而可以判断，网格4在一定程度上已经满足计算要求。

2.3 边界层网格的生成

边界层网格对定量的计算问题，如表面摩擦，升力系数阻力系数的求取影响比较大，对换热问题尤为显著。

尽管在工程精度允许的范围内，有时可以对复杂问题不画边界层，取而代之的是加密边界层网格，但由于本问题属于大雷诺数流动，且我们关注的就是阻力系数，因此有必要添加边界层网格，来提高计算精度。

Fluent中k-epsilon模型是针对充分发展的湍流模型，即高Re数的湍流模型。对近壁区内的流动，Re数比较低，湍流的发展并不充分，湍流的脉动影响不如分子黏性的影响大，这样，在这个区域内就不能直接使用该模型，通常采用壁面函数法或近壁模型法。

壁面函数法本身并不对贴壁层进行求解，而是用一组半经验公式来将壁面与充分发展的湍流区联系起来。其缺点在于，当流动情况偏离了壁面函数的理想条件时，如强体积力流动，近壁区三维性很强的流动问题，计算误差就会明显增大。其次壁面函数的缺点还在于：沿壁面法向细化网格时，会导致数值结果恶化。当y+小于15时，将会在壁面剪切力及热传递方面逐渐导致产生无界错误。

近壁模型法。修改湍流模型以使其能够求解近壁黏性影响区域，包括黏性底层。本文使用的方法即近壁模型法。近壁模型不需要使用壁面函数，但需要壁面网格的精细。 fluent中的k-omega湍流模型就是一种典型的近壁湍流模型。但要保证y+＜=1，以1为最佳。

根据给定的y+值，计算第一层网格与壁面的距离方法如下：

（1）计算雷诺数Re：

（2）计算表面摩擦系数：

（3）计算壁面切应力：

（4）计算摩擦速度：

（5）计算第一层网格高度：

继电器所使用的工质为25号变压器油，其密度ρ=895kg/m3，运动黏度ν=9.6×10-6m2/s。计算出动力黏度μ=νρ=8.592×10-3。计算中y+取1，特征速度U取2m/s。特征长度L取方形挡板的边长0.05m。代入式（1）～式（5）得第一层网格高度ΔS为8.1517×10-5m。

使用近壁模型法时，覆盖边界层的最小网格数量在 10层左右，最好能达到20层。还有一点需要注意的是，提高边界层分辨率，常常可以取得稳健的数值计算结果，因为只需要细化壁面法向方向网格。对于非结构网格，建议划分10～20层棱柱层网格以提高壁面边界层的预测精度。棱柱层厚度应当被设计为保证有15层或更多网格节点。另外，棱柱层大于边界层厚度是必要的，否则棱柱层会限制边界层的增长。这可以在获得计算结果后，通过查看边界层中心的最大湍流黏度，该值提供了边界层的厚度（最大值的两倍位置即边界层的边）。

平板湍流边界层厚度计算公式[7]为：

运动黏度ν=9.6×10-6，U=2m/s，L=0.05，代入式（6）得δ=3.4997452×10-3m。边界层网格厚度S应大于平板边界层取3.5×10-3。

假设网格按指数律分布，设总网格层数为n（=15），公比为r。

由式（7）可反解出公比r（=1.1374）。

至此，边界层参数首层网格高度h，公比r，层数n均得到，依据这些参数生成边界层棱柱网格。画好的网格如图5所示。

图5 计算网格剖面及边界层网格Fig.5 Compute Grid profile and boundary layer mesh

本文划分边界层后还有一个好处，由于本算例中继电器内部挡板是可以转动的，当挡板转动后，fluent的局部重构和弹性光顺算法动态生成的新网格，难以保证计算y+值的要求。所以定义边界层网格，并在 fluent将边界层指定为随同挡板一同运动的流体域。这就可以保证边界层不会因挡板的运动而破坏[8]。

2.4 边界条件

以入口管径80mm为特征长度，由式（1）可求出入口处雷诺数Re为10416.6，大于内部流动临界雷诺数Recr=2320，属于湍流，选用低雷诺数模型中的标准k-omega模型。

图中相应的边界均用文字加以说明。在 fluent中边界in表示进口。将其定义为 fluent中的速度入口条件velocity-inlet，取入口测试速度1.3m/s。Flunet对于亚声速流动，速度入口边界条件直接忽略压强的大小，故模拟中无压强条件，选用默认值即可。以下公式均来自ANSYS Fluent 14.5帮助手册。

（1）湍流强度I的计算：

（2）湍流尺度I的计算：

（3）湍动能k的计算：

（4）特性耗散率w的计算：

上式中为模型中一个经验常数其值近似为0.09。

依次求得式（9）～式（11）可得入口边界条件湍动能k=0.0152040。特性耗散率ω=40.2003909/m。

由于是不可压流体，出口速度和压强未知，出口边界“out”只能选用out flow条件。这实际上假定流场已充分发展，各个变量（压力除外）沿流动方向的梯度值为零。其他边界均定义为壁面。边界“wall”定义为默认的速度为零的壁面条件。

流场的计算采用SIMPLEC算法，梯度项、压力及动量方程均用二阶格式以提高精度。设置收敛残差为1×10-6。为提高计算精度，选用双精度求解器。

2.5 结果分析

首先看y+值分布。影响精度的只是挡板壁面的y+值，在这里首先观察挡板的y+值云图，如图6所示，左边是正面迎接冲击油流侧，雷诺数大，导致相应的y+值也大，右边是挡板背向油流侧，雷诺数小，相应的的y+值也小。

图6 y+值云图Fig.6 cloud of y+ value

为进一步衡量最大y+值是否满足低雷诺数模型k-omega模型的网格要求，即y+≤1。在挡板的正面取竖直线段up-down和水平线段left-right。分别绘制两条线段上y+值的散点图如图7所示。由图7可知，挡板上y+值除少数点大于1外，绝大多数均满足条件。所以挡板壁面处网格质量满足计算要求。

图7 线段left-right和up-down上y+值散点图Fig.7 y+ value on the scatterplot in line left-right and updown

在工程计算中对精度要求高时，网格和模型的选择尤为重要。本文参考文献[9]中给出第一层网格高度对敏感参数产生决定性的影响。本文结合理论计算来对比有无边界层网格的差异。

由于前边绘制的网格参考y+值是在特征速度为2m/s的条件进行的。由公式（1）～（5）可得，在一定范围内，y+值会随特征速度的增加而增加，为保证本计算y+值不会显著增高，满足y+＜=1的要求，我们分别在特征速度小于2m/s，取为1m/s、1.25m/s、1.5m/s、1.75m/s不同工况下，在 fluent中监控挡板表面摩擦系数。并将由式（2）计算的表面摩擦系数汇总于表2。

表2 表面摩擦因数计算结果

由图8数据整理可见，有边界层的数值计算结果和理论值吻合良好，没有边界层则使结果产生较大偏差。由此可见边界层网格对计算结果的影响很大。合理选择首层网格高度绘制的边界层网格，会大大提高计算的精度。

如图9所示，详细给出有边界层时，计算流场的细节。

图8 表面摩擦系数对比

图9 速度大小云图及流场的流线图Fig.9 Streamline velocity magnitude contours and streamline of flow field

3 结论

本文讨论了网格数量对CFD模拟精度的影响，用进出口压差作为对比参数，用倍增网格数的方法，对比不同网格密度对压差的影响，找到适合计算的最佳网格数。并严格按照计算模型的要求，绘制有助于提高计算精度的边界层网格，使壁面网格y+值满足模型准则要求。

在以表面摩擦系数理论计算值为参照，对比有无边界层网格时得到如下结论：边界层网格对计算结果的影响很大。合理选择首层网格高度绘制的边界层网格，会大大提高计算的精度。本文的意义还在于生成了分析计算所需的高质量的网格。为后续进一步计算（动网格计算）做好铺垫。

[1]康顺，石磊，戴丽萍，范忠瑶.CFD模拟的误差分析及网格收敛性研究[J].工程热物理学报，2010，12：2009-2013.

[2]刘重阳，于芳，徐让书.CFD计算网格误差分析的一个算例[J].沈阳航空工业学院学报，2006，04：21-24.

[3]杨永.基于网格计算的CFD模拟可信度分析[A].中国空气动力学会、中国宇航学会、中国航空学会、中国力学学会.计算流体力学研究进展——第十二届全国计算流体力学会议论文集[C].中国空气动力学会、中国宇航学会、中国航空学会、中国力学学会：2004：6.

[4]康顺，刘强，祁明旭.一个高压比离心叶轮的CFD结果确认[J].工程热物理学报，2005，03：400-404.

[5]晏丽琴，徐宏彤.基于CFD软件模拟精度与可信度影响因素分析[J].自动化与仪器仪表，2014，09：101-103.

[6]（澳）Jiyuan Tu Guan Heng Yeoh，（美）Chaoqun Liu 著，王晓冬译.计算流体力学——从实践中学习[M].第1版.沈阳：东北大学出版社，2009，10 ：147-148.

[7]陈卓如，金朝铭，王洪杰，王成敏.工程流体力学[M].第2版.北京：高等教育出版社，2004，1：360-361.

[8]隋洪涛，李鹏飞，马世虎，马富银，胡颖.精通CFD动网格工程仿真与案例实战[M].第1版.北京：人民邮电出版社，2013，05：92-94.

[9]纪兵兵，陈金瓶.ANSYS ICEM CFD网格划分技术实例详解[M].第1版.北京：中国水利水电出版社，2012，01：255-259.

[10]吴逸飞，邹正平.考虑来流边界条件不确定性的数值模拟方法[J].航空发动机，2015，02：35-39.

[11]康顺.计算域对CFD模拟结果的影响[J].工程热物理学报，2005，S1：57-60.

[12]刘智益.不确定性CFD模拟方法及其应用研究[D].华北电力大学，2014.

刘明亮（1990—），男，硕士，主要研究方向：流体机械内部流场分析，E-mail：1052746541@qq.com

李胜男（1971—），女，硕士，高级工程师，主要研究方向：电力系统继电保护等，E-mail：2450564588@qq.com

张思青（1957—），男，教授，E-mail：363685119@qq.com

Analysis of Mesh on the CFD Simulation Results

LIU Mingliang1，ZHANG Siqing1，LI Shengnan2
（1.Kunming University of Science and Technology，Kunming 650093，China；2.Yunnan Electric Power Research Institute，Kunming 650000，China）

In the CFD simulation process，the impact of highquality mesh on the calculation results is essential.This article describes in detail the process grid independence test.By Fluent double precision solver and k-omega model，the paper compares the calculation error of add and without boundary layer mesh，and then，Investigates the effect of the boundary layer mesh on calculation accuracy.

grid independent solution; grid independence test;k-omega models; boundary layer meshing

国家自然科学基金资助（编号：51369012）。