雷诺数为22000的二维方柱绕流仿真计算

1、15.1 RNG湍流模型在22000雷诺数二维方柱绕流中的应用【算例内容】：方柱尺寸为100mm×100mm，流场区域为1400mm×3500mm，方柱距离入口500mm。流体工质采用水，其密度为998.2千克每立方米，动力粘度为0.001003帕秒。由雷诺数为22000可以折算出入口流体的速度应为0.2210579米每秒。实验所用来流湍流强度为2%，取入口湍流长度尺度为0.007（入口湍流长度尺度取值参考文献十六P215）。【仿真目标】：对于雷诺数为22000的亚临界方柱绕流，合理选择湍流模型将是关键。对于分离流动，不适合采用标准k-湍流模型，在高雷诺数湍流模型族中，RN

2、G湍流模型是值得研究的选择。而SST k-模型、雷诺应力模型、分离涡模型和大涡模拟也是适合的，本节仅讨论RNG湍流模型。要求用两种近壁面处理方式处理方柱表面，第一种方式是采用尺度化壁面函数，第二种是采用加强型近壁面处理方法。计算完成后比较两种近壁面处理方法的结果。文献十七中给出了雷诺数为22000的方柱绕流的实验结果、大涡模拟结果和RNG计算结果。可以用来验证本节计算结果。关于雷诺数22000方柱绕流（关于壁面函数法、湍流粘度比的计算等概念可以参考文献十六）15.1.1 RNG模型与尺度化壁面函数在方柱绕流中的应用图15-1-1 分析流程图15-1-1给出了二维方柱绕流问题的分析流程，在Wor

3、kbench中建立模型，然后使用Mesh模块进行分网，最后提交到Fluent中进行计算。图15-1-2 流场模型按照算例内容建立流场模型如图15-1-2所示，然后进行分网，由于采用的是尺度化壁面函数法，所以近壁面网格没必要划分地太细。进入到Mesh模块后首先到mesh模型树下将Physics Preference改为CFD与Fluent，然后设置Sizing下的网格相关度、网格光顺等参数。图15-1-3 网格划分结果图15-1-3给出了网格划分结果，在划分网格时，约束整个面的网格尺寸为10mm，并且定义入口、出口、上下边界、方柱表面的名称集（Named Selections），以方便在Flue

4、nt中定义边界条件。图15-1-4 局部网格示意图15-1-4给出了局部网格示意图，从图中可以看出，在方柱表面的网格高度是较高的，使用尺度化近壁面函数法时并没有划分很细的边界层网格。最后共得到48900个网格，网格最大长宽比为1.0069，满足仿真计算的网格要求。图15-1-5 采用4核并行计算在打开Fluent时，按照图15-1-5所示，采用单精度、4核并行计算。双精度一般用在以下几个方面：（1）存在明显的细长比很大的几何特征，如细长的管道；（2）流体是几个区域，区域之间用小尺寸的管道连接，其中一个区域的流体压力大大高于整个流域的平均水平；（3）对于存在较高导热率的共轭传热问题；（4）对于流

5、场网格长宽比较大的问题。在本算例中，虽然不存在以上问题，但仍采用双精度4核并行计算求解。在进行计算时，先进行稳态计算，然后将稳态计算的结果作为瞬态计算的初场，再进行瞬态计算。求解器采用压力基求解器，并进行湍流模型的设置。图15-1-6 湍流模型的选择如图15-1-6所示，选择RNG湍流模型与尺度化壁面函数处理方法，并且不启用Differential Viscosity Model选项，在ANSYS Fluent 14.5的帮助文件中对该模式有以下论述：The RNG turbulence model in ANSYS FLUENT has an option of using a d

6、ifferential formula for the effective viscosity to account for the low-Reynolds-number effects.可见在考虑到边界层的低雷诺数影响时应该打开该选项，另外，有人在做空气的涡街仿真时遇到过不勾选该选项就得不到涡街的情况。（但是在实际计算中发现，启用该选项后，计算容易变得不稳定，出现了数值波动甚至是不收敛。因此，在本算例中，不启用该选项。流线曲率修正是因为使用各向同性涡粘假设的湍流模型对流线曲率不敏感，因此在流线弯曲时需要进行流线曲率修正，启用该选项。选择Fluent中的water-liquid作为流体材料，

7、然后进行边界条件的设置：设置方柱表面为无滑移壁面条件wall；流场上下边界为对称边界条件symmetry（在数学上即自由滑移边界条件）；当流场发展到出口位置时，认为流动已经充分发展，因此设置出口为自由流outflow；入口在用速度入口边界条件，入口速度、湍流强度、湍流长度尺寸均采用算例中所给设定。在设置参考值时首先选定compute from->流场名，这样就正确的给出了参考值中的物理参数（密度、粘度等），然后指定Area的值为0.1。关于参考值的设定是非常有讲究的，根据文献十八P173的升阻力系数定义式：Cd=Fd12v2SCl=Fl12v2S以上两式中，Cd代表阻力系数，Cl代表升力

8、系数，Fd为钝体所受到的阻力，Fl为钝体所受到的升力，为流体的密度，v为前方来流的速度，S代表钝体的投影面积，在二维问题中S的数值与弦长相等（弦长乘以一，二维问题中展向长度为单位长度一）。有关翼型的概念：翼型：飞机机翼具有独特的剖面，其横断面（横向剖面）的形状称为翼型，称为翼型前缘：翼型最前面的一点。 后缘：翼型最后面的一点。 翼弦：前缘与后缘的连线。 弦长：前后缘的距离称为弦长。如果机翼平面形状不是长方形，一般在参数计算时采用制造商指定位置的弦长或平均弦长攻角：也称迎角。对于翼形来说，攻角定义为翼弦与来流速度之间的夹角，抬头为正，低头为负，常用符号表示。其示意图见图15-1-7：图15-1-

9、7 机翼攻角示意图升力系数：物体（如飞机、导弹）所受到的升力与气流动压和参考面积的乘积之比，是一个无量纲量。阻力系数：对于飞行器来说，阻力系数定义为物体（如飞机、导弹）所受到的阻力与气流动压和参考面积之比，是一个无量纲量。关于升力系数与阻力系数的计算参考图15-1-8：图15-1-8 计算升阻力系数所用的参考系需要注意的是，在计算升力与阻力时，升力方向垂直于来流方向（不是垂直于机翼！），阻力方向平行于来流方向（不是平行于机翼！）。仅当攻角为0°或180°时，升力/阻力才垂直/平行于来流方向。在Fluent二维计算中，机翼长度方向（展向）垂直于几何平面，软件在处理时机翼长度置

10、为1（在流体力学中处理二维问题，展向长度全部取单位长度）。其余参考值需要参考升力系数公式进行计算，并且一定要注意升力/阻力的方向。对于本二维算例，S的取值与弦长数值相等，因此本算例S取值为0.1，即设定参考值时，与尺寸相关的进设置Area为0.1平方米即可，其余的保持默认不变；密度等材料的参考值已经通过Compute from功能实现了。之后设置求解器选项，选择耦合求解器，对流项的处理采用二阶迎风格式（对于之前划分的网格，也可采用更高阶的QUICK格式）。保持默认的收敛残差判定标准0.001不变，监测了升阻力系数随迭代步数的变化曲线。以入口条件初始化整个流场区域，设定最大迭代步数为1000步。

11、图15-1-9 稳态计算的残差曲线图15-1-9给出了稳态计算残差曲线，从曲线中可以看出，残差曲线基本上随着计算会继续趋于收敛，其中连续性方程的残差最大。在计算此类绕流问题的实际应用过程中，还会出现当稳态计算的残差位于高位稳定的情况，出现这种情况时，应当考虑的一种可能存在的情况就是物理模型中存在瞬态变化的物理效应，导致稳态迭代时前后迭代的范数残差不能够收敛而位于高位稳定，这种现象在用稳态雷诺平均方法求解瞬态问题时经常出现，此时可以将稳态计算结果作为瞬态计算的初场，继续计算。所以单看残差收敛情况不一定有太大意义，要结合实际情况来看。另外发现当启用Differential Viscosity Mo

12、del选项后，本来收敛至非常小的残差容易转变为高位稳定甚至发散。本例中升力系数与阻力系数的迭代历史曲线如图15-1-10所示：图15-1-10 方柱绕流的稳态升阻力迭代历史从图15-1-10中可以看出，升阻力系数在稳态计算中最终归于稳定，这与残差收敛至极小值是相互呼应的。而在有些绕流问题的稳态计算中，也会出现瞬态特征。如图15-1-11所示： 图15-1-11 某圆柱绕流的稳态计算图15-1-11中给出了某圆柱绕流的稳态计算结果，在稳态计算结果中出现了瞬态特征。在稳态计算中出现瞬态特征的现象，在方柱绕流与圆柱绕流的现象中都可能出现，当出现这一现象时，残差曲线将会居于高位，并且升力与阻力系数的值

13、不一定准确，比如上图中升力系数平均值部位0。还有一种现象是升阻力系数趋于平稳，但残差曲线局于高位，造成这种现象的一种原因是启用了Differential Viscosity Model选项，这也说明了Differential Viscosity Model选项对计算所造成的不稳定性。根据文献十八P119的关于斯特劳哈尔数的公式：St=fDv其中，St为斯特劳哈尔数，f为漩涡脱落频率，v为来流速度。根据文献十七中的实验结果（St=0.13），可以反推出漩涡脱落频率应为3.5s，所以在选择瞬态步长时就有了选择的依据，瞬态求解采用二阶隐式求解方式，针对方柱绕流问题，隐式方法对时间步长的要求不高，很多

14、人取漩涡脱落周期的1/50作为瞬态计算的时间步长，本文取为0.02s，可以充分地描述出漩涡脱落的瞬态细节。对于收敛残差的取值，很多人是取稳态残差稳定值的1/50，本文取为1/100，即0.00001。指定每一时间步最多计算25个子步，共计算5000个时间步，共计100s。图15-1-12 方柱绕流流场的静压云图与动压运图根据图15-1-12中可以看出，在方柱后方出现了明显的周期性漩涡脱落的现象。漩涡脱落的现象在总压云图与速度云图中也可以很明显地看出来。文献十七中认为RNG方法在预测平均阻力系数和回流区再附着距离时与实验及大涡结果比较接近，但是在预测升阻力系数的均方根值时仍具有一定的差距。采用雷

15、诺平均的RNG模型，从模型精度上来说要比大涡模拟至少低一个数量级。图15-1-13 方柱绕流流场的总压云图与速度云图图15-1-14给出了方柱升阻力系数的历史曲线。图15-1-14 方柱绕流的升力与阻力系数的历史曲线从图中可以明显看出，方柱绕流的升阻力系数随着时间的推移呈周期性的波动。经过对数据的统计，可以知道，方柱阻力系数的平均值为2.3（取80-90s之间的每个时间步的阻力系数做平均），实验值为2.1，误差为9.5%，计算出的漩涡脱落周期为3s，由此所得的斯特劳哈尔数为0.15，实验值为0.13，误差为15.4%。图15-1-15 方柱表面的Y+值分布范围仿真后可以知道，近壁面Y+值的分布

16、范围为7.07-68.6，在Fluent的官方文献中推荐壁面函数法的Y+值适用范围为30-500，对于尺度化壁面函数，由于在近壁面使用了预防Y+值过小导致结果恶化的技术，其Y+值范围下限可以适当放宽，但仍应尽量取在11.25以上。过大或过小的Y+值都可能引起结果的恶化。15.1.2 RNG模型与加强型近壁面处理在方柱绕流中的应用图15-1-16 加强型近壁面处理所用到的流场网格返回到MESH模块中，将方柱边的尺寸约束为1mm，在方柱的边上添加边界层网格控制，指定第一层网格高度为0.25mm，网格增长率为1.05，共划分了15层边界层网格，得到了135500个节点，共计68053个单元。所有单元

17、的最大长宽比为4.53，完全满足加强型近壁面处理的网格计算要求。图15-1-17 加强型近壁面处理后的方柱升阻力系数历史历程通过图15-1-17可以看出，使用了加强型近壁面处理方法之后，所得到的升阻力系数波动更加稳定直观。将升阻力系数的保存文件使用Matlab读入并处理后，可以知道，加强型近壁面处理计算出的波动稳定后的阻力系数平均值为2.15，与实验值2.10的误差为2.38%，准确度明显高于尺度化壁面函数法；计算所得的斯特劳哈尔数为0.14，与实验值0.13的误差为7.69%，其准确度也明显高于壁面函数法。计算所得的阻力系数波动的均方根值为0.125，而W.Rodi等人通过大涡模拟得到的阻力

18、系数波动的均方根值为0.14，误差为-11.4%（文献十七使用RNG模型计算的结果为0.11）；计算所得的升力系数波动的均方根值为1.41，W.Rodi使用大涡计算的结果为1.01，误差为39.6%（文献十七使用RNG模型计算的结果为1.26）。通过以上结果可以看出，RNG模型与加强型近壁面处理方法的结合使用，能够较好地预测阻力系数平均值、漩涡的脱落周期与阻力系数波动的均方根值，但是在预测升力系数波动的均方根值时，还有较大的误差，与大涡结果存在一定的差距。（注意，升阻力系数波动的结果均是与大涡模拟比较，在查阅资料时并未发现相关实验数据，因此在评判RNG计算所得的升阻力系数波动时，均是与大涡模拟

19、结果相比较）同时必须检查Y+值是否处于合理的范围之内，通过Contour中的Turbulence->Y Plus来显示方柱表面上的Y+值分布范围，经查看后可以知道，Y+值分布在0.197-0.547之间，稳态计算结果中的Y+值比瞬态计算的Y+值略小，但是前后两次计算的Y+对于RNG与加强型近壁面处理来说都是可以接受的。对于静压、动压、总压、速度云图如下图15-1-18所示：图15-1-18（a） 静压与动压云图图15-1-18（b） 总压与速度云图附录：简单Fluent升力、阻力、力矩系数后处理程序%简单Fluent升力、阻力、力矩系数后处理程序 徐太栋 2013clcclearformat longfid=fopen('S