Cd计算不准的问题-LES-Re3900-三维圆柱绕流-pisoFoam


  • OpenFOAM副教授

    各位OFer:

    最近我在做三维圆柱绕流问题,但是Cd却无法和文献出现的数据相匹配,问题不知道出在什么地方。
    下图中 蓝色曲线是drag系数历史,20个周期的平均值为1.3左右,Strouhal 0.20和文献非常接近。但是Cd却得到一个偏高的数值。0_1482344985625_Result - Copy.png

    网格是用ICEM按照Youtube上的视频一步步做的。首先生成的是二维网格,然后在span wise 的方向上伸长了32个节点,长度在 pi*D/2。 我用2维网格在ANSYS/FLUENT 17中做过测试,Re在3900的情况下,能满足Y+小于1。网格形态和文献中出现的矩形情况类似。 32节点的三维扩展,得到了的总网格数在5百万左右。 之后,我又做过补充测试64个span wise 的伸长,但是结果没有太大差别。
    1_1482345271911_Y+.JPG 0_1482345271911_Grid2.JPG

    边界条件。 速度INLET,零压力OUTLET, noSlip 圆柱表面。 上下面是墙函数并且使用零剪应力条件,前后面是对称边界条件。之后,我还尝试了cyclic边界对于上,下,前和后面,结果还是没有太大差别。

    时间步长取的是10^-4秒,求解器pisoFoam, SpalartAllmarasDDES LES 湍流模型,OpenFOAM v4.1 版本,硬件和操作系统是:24 cores/node, Memory per node 32 G, infiniband, AMD @2.1GHz CPU, Centos 6.8. 大部分是用144核并行计算。

    0,constant和control文件如下:
    10_1482345647509_U.txt 9_1482345647509_turbulenceProperties.txt 8_1482345647509_transportProperties.txt 7_1482345647509_p.txt 6_1482345647509_nuTilda.txt 5_1482345647509_nut.txt 4_1482345647509_k.txt 3_1482345647509_fvSolution.txt 2_1482345647509_fvSchemes.txt 1_1482345647509_decomposeParDict.txt 0_1482345647499_controlDict.txt

    大家有什么样的建议吗?



  • 你的进口是否已经很长了?SpalartAllmarasDDES LES换过别的么?


  • OpenFOAM副教授

    inlet 到圆柱中心距离是10倍的圆柱直径,outlet到圆柱中心是20D, 上下两个面距离是20D。

    之前在FLUENT上用同一套网格做过3900的例子,Cd均值为1.1,WALE是LES的子网格。

    别的子网格还没有尝试,我会尝试用别的子网格。

    现在500万的网格正在尝试用更严格的收敛准则,10的负7次方对于所有变量。



  • @random_ran 圆柱长度用piD应该好点,你看看你算出入口的p是不是大?如果大就是入口和出口之间存在的压差对结果造成了影响,如果大就要对入口到圆柱之间的距离加大,或者适当调整一下网格。还有你看看你的courant number,我感觉这个不要太大,同样网格下,不同的courant number结果也会有不同。


  • OpenFOAM副教授

    @sorrby

    piD 在昨天跑完10个周期之后,发现和piD/2,没有太大差别,1.3附近,稍小一点。

    Courant number在整个计算中保持在0.2以下。

    对于压差问题:

    1. 出口边界条件是 fixValue, 数值是0;
      2.入口的压力边界条件设置是 zeroGradient;

    流场还在重构中,传回我电脑还需要一点时间。等结果出来之后,我会回复你的。

    谢谢你的回复。



  • @random_ran 不行,你就在试一试换一种网格画法,调整一下你的网格。


  • OpenFOAM副教授

    @sorrby

    对于这个网格,个人感觉是能接受的。主要原因是用同样的网格用FLUENT跑出了1.1的Cd。

    到目前为止,我其实还有很多地方不太明白。特别是边界条件的设置,离散格式的选取。

    对于边界条件:
    OF非常不友好。要在0文件中设置5个子文件。每一个子文件中k,p,U,nut 以及nuTilda 对每一个边界都要从头到尾设置。 这套边界条件的设置,特别是k,nut和nuTilda是基于教程里的摩托车算例。但是对于他们的物理意义理解就不如p和U那么明朗。比如k文件中的0.24取值:在inlet是fixedValue,在outlet是inletOut边界条件。 还比如nut和nuTilda文件的设置。我不清楚这些参数会对Cd有什么样的影响。

    另外对于cyclic 边界条件,我尝试过的是上,下,前 和后面,四个边界条件同时改成cylic,会不会是只用使用前后两个面cyclic,另外两个面用wall+zeroGradient?

    对于时间的离散,我用的是一阶精度的backEuler,这个对比FLUENT有所不一样,FLUENT会主动提示你时间离散选用二阶以上精度。问题可能出在这么?

    另外还有一个比较极端的想法,既然FLUENT能算出1.1,OF用同样的设置不行吗? 问题在于FLUENT的LES子网格我使用的WALE,OF的子网格有很多,难道说应该把主要精力放在LES的子网格上吗?


  • OpenFOAM教授

    backward 是二阶精度,Euler 是一阶精度。LES 和 DES 建议采用二阶精度的 backward。

    backwardEuler 是啥?

    Re=3900 的圆柱阻力系数的实验值应该是 0.99±0.05,1.1 差的还是有点远。


  • OpenFOAM副教授

    @wwzhao

    不好意思,是我犯错了。 fvScheme中用的是二阶精度的backward。


  • OpenFOAM副教授

    @sorrby

    的确,在从o-shape形状区域到wake区域的转换有一点瑕疵。转换比太大?我会重画一下网格再试试。


  • OpenFOAM副教授

    首先,祝大家新年快乐! 对一个月前的问题继续进行补充。

    用了 dynamicKEq 的子网格,得到的阻力和升力系数如下:

    0_1485887503622_Cd-dynamicKEqn.png
    0_1485887533502_Cl-dynamickEqn.png

    计算还在继续,虽然现在只算了几个周期。Cd的均值在0.99左右,这个值和文献报道的数据非常接近了,但问题也很明显:居然出现了负的阻力系数。之前的模拟虽然Cd稍微偏高,但是阻力系数没有出现负值,甚至都没有下过1.00。网格是用的同一套,y+在1以下,Courant 最大值在整个过程中在0.2以下。

    Lysenko, D. A. et al.(2012)同样用OpenFOAM进行了LES模拟三维圆柱绕流,其中他对比了两个子网格的不同,传统的SMAG以及dynamic k-equation(TKE)。论文中的Cd,Cl 历史如下:
    0_1485884234007_reference.png

    他们的TKE模型计算得到的Cd的r.m.s 值很小,和我的试算反差很大。查看fvSolution文件后,我怀疑是残差标准的设置可能会影响数值计算结果,所以会想办法提高标准后看结果怎么样。

    另外一个比较头疼的问题是尾流的速度曲线:文献中出现的这根曲线,横轴是计算域中轴线(原点在圆柱中心);纵轴是进口速度方向上的速度大小,通常是非常平滑的,如下图:
    0_1485889453783_streamwiseVel.png
    而我的计算结果,请注意是黑色虚线那条:
    0_1485889536137_recirculationLength.png

    虽然大体上类似,但是波动却非常明显,不知道这是什么原因造成的。我的采样点是从Cl的历史中,发现振荡稳定后,取了10个周期。 而上篇文献则取了150多个周期,这个问题的来源吗?或许大家有更好的见解?


  • OpenFOAM副教授

    换了计算域的形态,居然就能吻合得很好了。
    Cd 均值 1.03; St 0.21; Cl_rms 0.20; 30个周期96核,只用了12个小时;
    可是为什么之前矩形计算域,算出偏高的数值?
    难道雷诺数一样,不同的流体性质,也会对计算结果产生影响?

    2_1487101846107_compDomain.png

    3_1487101846107_Cylinder.png

    4_1487101846107_Mean stream-wise velocity.png

    1_1487101846107_Cl-dynamickEqn.png

    0_1487101846107_Cd-kEqn.png



  • @random_ran 你这次只改了流场吗?你之前矩形的流畅有没有把过度的地方改一改,看看结果有影响吗?相同雷诺数,不同的流畅性质应该不会产生影响。


  • OpenFOAM副教授

    @sorrby 改了很多地方,比如fvm 的scheme, 时间步长,收敛标准。

    对于之前那套网格,实在是没有精力了。一个二维的网格就16W的cell了,已经做得很密了。 反倒是圆形的这个计算域,只用到4W的cell就把很多参数算得很准了。要是有时间,再试试那个矩形的计算网格。

    另外,我的意思是这样的,对于相同的雷诺数,不同的流体性质,会不会对一些测量参数产生一定范围内的影响?比如10%以内。



  • @random_ran 做科研的态度:laughing:


  • OpenFOAM副教授

    @李东岳

    多谢前辈鼓励!

    Proud to be a FOAMer.

    :happy:


  • OpenFOAM讲师

    你好,我做了一个类似的算例,不过是二维的试算,我查看了涡量,不知为何非常大,请问你有没有类似的问题呢?


  • OpenFOAM讲师

    另外还有,请问你的时间平均怎么做的呢?比如第一个图蓝色的线Cd,是把每个周期内t~t+delta t的Cd做平均画出的图,还是直接吧0~t时间段的做平均?


  • OpenFOAM副教授

    @yhdthu

    Hi yhdthu:

    我的这个case没有重点关注涡量的绝对大小。主要考察的是平均Cd, Strouhal 数,recirculation的长度。 voriticiy 只是用三维等值面做了一下wake的结构。 不知道你的case 重点考察的对象是什么? 如果你能提供的一下你的涡量考察的座标,我可以把我的数据提供给你。

    对于蓝色的线的Cd 我是把每个一个时间步长的Cd都输出来。这个Cd是平均的整个Cylinder的表面。 然后把每个时间点的Cd 对时间用matlab plot一下就是那张图了。


  • OpenFOAM讲师

    @random_ran
    好的,谢谢你
    第一个应该是用的DES模型计算的吧,我看帖子最上面的Cd随时间变化幅度很小,但是后面哪个用了动力模式做的Cd图波动很大,虽然均值基本上是和实验对得上的

    我的目的是想做空化模拟,但事前要先将湍流结构捕捉的较好才可以,目前二维试算用的是动力模式,Cd波动也很大,我不是每个时间步输出的结果,平均起来效果不是很好,可能是样本点太少了

    可否将Cp曲线也po出来呢?还有其时间平均的定义,希望能做个参考,非常感谢

    BTW,如果可以的话,加个微信好友,可以随时交流~358253794


  • OpenFOAM副教授

    最开始确实是DES。 dynamicKEqution 的一大特点就是波动性,我第一次看到这个湍流模型运行出来的Cd历史的时候也被惊讶到了。可以参考:

    Lysenko, Dmitry A., Ivar S. Ertesvåg, and Kjell Erik Rian. “Large-Eddy Simulation of the Flow over a Circular Cylinder at Reynolds Number 3900 Using the Openfoam Toolbox.” Flow, Turbulence and Combustion 89.4 (2012): 491–518. Web.

    如果湍流结果是一个前置条件的话,特别是你的雷诺数上千以后,2D的模拟总会计算出一个较高的Cd结果。 我不知道你是怎么计算Cd,O.F. 自己有专门计算Cd的object:

    FoamFile
    {
        version     2.0;
        format      ascii;
        class       dictionary;
        object      forceCoeffsIncompressible;
    }
    

    你可以按自己的需求调整输出频率。另外,文献中对采样点的多少也是各有各的看法。有的人(Franke 2002)采样150个周期(shedding period), 有的人(Kravchenko and Moin 2000)说采样10几个周期的结果和自己的实验结果很吻合,而更多的周期则会出现偏差。个人建议30~50个周期。

    我这里有最开始那个DES的Cp图,传统的plot方法。后面的因为考察的方向不同,没有用传统的方式画,对你也没有参考价值。Lysenko 2012 年那篇里有Cp的plot,你可以参考他的。

    Cp的做法,我是用ParaView 后处理得到的。先在OF设置好:在一个周期内,输出,比如20个点 (包含p,U场)。计算完成后在圆柱的中截面截取出压力值用来计算Cp (这个Cp是截面Cp值),这和Cd (整个圆柱)的算法是有区别的。得到了一个个的输出的时间点的Cp值之后,再对所有的采样的周期,比如30个周期,求出平均值。

    ps 论坛的回复没有办法插入图片,如果你有需要我那个Cp,可以发我的邮箱.

    个人习惯:不用微信以及一切即时通信软件。这是并行计算的哲学告诉我的:尽量减少cpu之间的通信以提升效率。 :joking:


  • OpenFOAM讲师

    @random_ran 好的,方法基本一致,我自己再整理一下,我觉得我的问题应该是后处理的问题,有问题再交流,祝好


  • OpenFOAM教授

    @yhdthu 可以做fieldAverage,然后用生成的pMean来计算Cp。


  • OpenFOAM讲师

    @wwzhao 是的,我就是这么做的


  • OpenFOAM讲师

    还有个问题,请问LES网格大小是怎么确定的呢?我是根据kolmogorov理论确定的delta,即:

    L0/6 < delta < 60L0Re^(-3/4)

    (L0是积分尺度)

    还有壁面附近是怎么处理的呢?

    依据此做的网格,感觉网格有点粗,不知你是怎么确定的呢?


  • OpenFOAM副教授

    @yhdthu

    近壁的地方是最难处理的,因为在这些地方没有大涡。LES也可以用wall function来避免太多的计算量,但是wall-resolved的LES的网格和DNS非常接近。

    通常对壁面网格的控制是用无量纲墙距离 (yplus/y+)来标定。这个量的确定和雷诺数密切相关。对于wall-resolved的LES,y+<1是文献中最常提到的。对于用wall function的例子,y+可以放到30到200之间用来避免buffer layer,不同的wall function 还是有一些不同的要求。

    对于stream-wise, span-wise,wall-resolved 的LES还有更高的要求。Georgiadis (2010)给出的建议是:

    50< Δx+<150; Δy+<1; 15< Δz+<40

    POINTWISE, NASA都给出了依据y+计算第一层网格高度的公式。

    可以看看这个帖子


  • OpenFOAM讲师

    对于方形计算域,我猜是不是因为网格过渡均匀性做的不好导致的Cd偏大?



  • 方形计算区,o型网格覆盖3D~4D范围,再加上尾迹区加密,结果应该是一样的。主要原因,猜测是非正交网格对涡解析不够,造成Cd或Cl差异。至于涡对Cd和Cl影响,可以参考吴介之老师的涡动力学。


  • OpenFOAM副教授

    当时的计算结果高了大概10%~20%,当时试了很多:
     Fixed grid jump transition area
     Use LUST scheme for the convective terms
     Implement LES turbulent models such as TKE and Smagorinsky model
     Use nutLowReWallFunction instead of nutUSpaldingWallFunction
     Increase computational domain
     K and nut parameters
     Boundary conditions: two cyclic with two slip wall conditions

    最后的结果发现换网格的效果是最好的。圆形计算域的正交性更好一些。而且这个计算域用自动生成都可以做到很好的质量,何乐而不为呢?



  • @Junhua-PAN 求书名,网上没搜到



  • @random_ran 想请教一个问题,如果用symmetry网格的话,对计算结果影响大吗?


  • OpenFOAM副教授

    如果单纯从均值Cd的大小上看,这两个边界条件,对于我所用的那个网格并没有太大影响。

    我最开始的时候也是用过symmetry的边界条件,但是我觉得slip-wall 的边界条件可能更适合描述我当初对这个问题的理解。

    后来我发现这个研究方向上的文献中,绝大多数人在span-wise都是用cyclic的边界条件。用cyclic边界的条件的优势之一就是能把模拟的有限长度的圆柱,当成无限长来出来。换句话说就是end-effect上有很强的优势。这个即便是做风洞实验都无法克服的问题,因为风洞中都要用end-plate来矫正端部对流场的影响。

    Yeo, DongHun, and Nicholas P. Jones. “Investigation on 3-D Characteristics of Flow around a Yawed and Inclined Circular Cylinder.” Journal of Wind Engineering and Industrial Aerodynamics 96.10–11 (2008): 1947–1960. Web.

    这篇文章中3.1节讨论了 slip-wall 和 cyclic wall 对涡结构的影响。作者发现在cyclic边界条件受边界的影响要小,因此cyclic边界条件能更好模拟三维圆柱绕流这一经典问题。



  • @random_ran :kiss: :kiss:



  • @yuan_neu Vortical Flows 第9章



  • @random_ranCd计算不准的问题-LES-Re3900-三维圆柱绕流-pisoFoam 中说:

    时的计算结果高了大概10%~20%,当时试了很多:
     Fixed grid jump transition area
     Use LUST scheme for the convective terms
     Implement LES turbulent models such as TKE and Smagorinsky model
     Use nutLowReWallFunction instead of nutUSpaldingWallFunction
     Increase computational domain
     K and nut parameters
     Boundary conditions: two cyclic with two slip wall conditions
    最后的结果发现换网格的效果是最好的。圆形计算域的正交性更好一些。而且这个计算域用自动生成都可以做到很好的质量,何

    你对比的是啥文献,用的什么边界条件呀?

    我看有个早期文献是用FVM加变换算的无穷大边界的,和一般的有限大小边界还不一样。

    如果你考虑正交性,你在fvScheme中加修正了么?梯度用的啥?OF对非正交和skew都有相应的修正方法,可以提高一点儿精度。我写过一个小帖子OF修正小结,你可以尝试搞搞边界条件类型和fvScheme中的修正。


  • OpenFOAM副教授

    主要参考的是
    Lysenko, Dmitry A., Ivar S. Ertesvåg, and Kjell Erik Rian. “Large-Eddy Simulation of the Flow over a Circular Cylinder at Reynolds Number 3900 Using the Openfoam Toolbox.” Flow, Turbulence and Combustion 89.4 (2012): 491–518. Web.

    Lysenko (2012) 用的是:

    • inlet: 层流
    • outlet: wave-transmissive conditions (不太明白,我用的是 pressure:fixValue=0)
    • span-wise: periodicity
    • cylinder: no-slip

    能提供一下那篇早期文献吗?

    当时没有考虑正交性问题, 关于梯度的离散是这样设置的:

    gradSchemes
    {
        default         Gauss linear;
        grad(nuTilda)   cellLimited Gauss linear 1;
        grad(U)         cellLimited Gauss linear 1;
    }
    

    这个贴子当时是在算Re3900的时候遇到Cd的过高估计,最终发现换一套网格和用dynamic kEquation 就能把Cd预测的很准了,所以那段研究就没有再进行下去了。

    有时间我会再把以前的那套网格用你提供的方案再算算,看看结果会怎么样。



  • @random_ran
    waveTransmissive类似特征边界条件,反正Wikki是这么说的。可能更加有物理意义吧。


  • OpenFOAM讲师

    你好,我用DDES-SA模型做出来的模拟图如下:
    0_1494666752922_WechatIMG8553.jpeg

    这个对的不是很好,另外流向平均速度就更离谱了,选取了流向三个位置,
    x / D = 1.06;x / D = 1.54; x / D = 2.02,如下

    0_1494666863390_WechatIMG8554.jpeg

    我对比了两个文献,其中一篇是你推荐的,如下:
    0_1494666927269_WechatIMG8555.jpeg

    0_1494667023928_WechatIMG8556.jpeg

    我不明白,都是无量纲分布,怎么差别这么大?我的计算问题会出在哪呢?谢谢



  • @yhdthu

    模型与模型之间可能是不一样的,参考NASA TMR 关于SA模型的实现,https://turbmodels.larc.nasa.gov/spalart.html#sa,SA模型有10种左右变种。

    OF实现的是SA-fv3.


  • OpenFOAM讲师

    @程迪 你好,我用的模型是DES类模型,在我取点处已经切换成LES模型计算了,另外无量纲分布怎么会随着模型变化呢?把它变成有量纲形式就可以看出来不可能,分布是一定的,可以参考实验


  • OpenFOAM讲师

    @yhdthu 经过查阅文献,确认了文献纵坐标是废的。。


  • OpenFOAM副教授

    @yhdthu

    我觉得你第一个图做得已经很棒了!比对的是两组实验数据,你可以再加上数值模拟方面的数据就更好了。

    如果你想再提高精度,更精密的网格肯定逃不了。另外湍流的模型的选择对结果的影响很大。还有就对流项的离散选择。关键是看你模拟的目的是什么。

    Ps: 文献纵座标是废的是什么情况?能详细说说吗?
    Ps2: 网格数量是多少?求解核心数?计算时间?边界条件如何?不知能否分享一下你的经验?


  • OpenFOAM讲师

    @random_ran 你好,我的意思是纵坐标其实没有实际意义,可以想象,在y/D>=3时,Ux=Uinf,所以文献给出的其实是示意图,是为了方便说事。
    另外,我用的就是方形计算域,大概是220W,壁面附近特殊处理了一下,分成三个增长率生成网格,过渡区网格注意了一下。基于升力的St大概是0.202,Cd约为1.02,下次在电脑旁我会把其他参数加上。
    对流项离散LUST,就是那个75%CD,25%UD
    DES类湍流模型,DDES_SA


  • OpenFOAM讲师

    还有个问题,我想做出平均分离角度,使用的是哪个类呢?我看有个wallshearshtress,但是这个是湍流雷诺应力R与壁面法向量n的内积,应该不是,请问你是怎么处理的呢?


  • OpenFOAM副教授

    wallShearStress可以在计算完成后处理,O.F. (v4.1)是有这个功能的,但是个人觉得这个没法找到分离点的座标,只是给出了应力的最大和最小值。

    simpleFoam -postProcess -func wallShearStress
    

    在分离的地方,圆柱表面的剪切应力应该很接近零。不过一个点的剪切应力应该有三个分量,是考虑某个分量上的最小值还是合成的最小值?我也在想这个问题。

    不过wallShearStress是和壁面速度沿壁面法向的梯度,所以用这个梯度也可以来求分离点的位置。结合ParaView是很容易实现的。

    我建议还是用ParaView 后处理。

    Filters-> Alphabetic->Gradient of Unstructured Dataset 对速度场进行处理。

    然后用 Plot on intersection curve 输出圆柱表面的梯度分布,再输出所有时间步长的数据,就可以求出时间上的均值。如果在多切几个不同的界面,就可是实现空间上的均值。

    希望能帮助到你。


  • OpenFOAM教授

    @yhdthu 分离角可通过gradU计算得到。


  • OpenFOAM讲师

    @wwzhao 你好,可以说说具体内容不?谢谢


  • OpenFOAM教授

    @yhdthu 上面说错了,应该是用wallGradU计算得到壁面处的速度梯度,du/dn=0 的点(n为圆柱表面法向量)即为分离角所处位置。


  • OpenFOAM讲师

    @wwzhao 你好,我试过这个,在paraview里显示的是magnitude,并没有0的点,请问怎么检测呢?

    0_1495124636735_WechatIMG8730.jpeg


  • OpenFOAM教授

    @yhdthu

    1. 你需要的是平均分离点,因此需要修改wallGradU的源码,让它读Umean,而不是U。

    2. 你需要做一个切片,对切片上的wallGradU进行分析。


Log in to reply
 

CFD中文网 2016 - 2020 | 京ICP备15017992号-2