Cd计算不准的问题-LES-Re3900-三维圆柱绕流-pisoFoam
-
各位OFer:
最近我在做三维圆柱绕流问题,但是Cd却无法和文献出现的数据相匹配,问题不知道出在什么地方。
下图中 蓝色曲线是drag系数历史,20个周期的平均值为1.3左右,Strouhal 0.20和文献非常接近。但是Cd却得到一个偏高的数值。网格是用ICEM按照Youtube上的视频一步步做的。首先生成的是二维网格,然后在span wise 的方向上伸长了32个节点,长度在 pi*D/2。 我用2维网格在ANSYS/FLUENT 17中做过测试,Re在3900的情况下,能满足Y+小于1。网格形态和文献中出现的矩形情况类似。 32节点的三维扩展,得到了的总网格数在5百万左右。 之后,我又做过补充测试64个span wise 的伸长,但是结果没有太大差别。
边界条件。 速度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换过别的么?
-
inlet 到圆柱中心距离是10倍的圆柱直径,outlet到圆柱中心是20D, 上下两个面距离是20D。
之前在FLUENT上用同一套网格做过3900的例子,Cd均值为1.1,WALE是LES的子网格。
别的子网格还没有尝试,我会尝试用别的子网格。
现在500万的网格正在尝试用更严格的收敛准则,10的负7次方对于所有变量。
-
@random_ran 圆柱长度用piD应该好点,你看看你算出入口的p是不是大?如果大就是入口和出口之间存在的压差对结果造成了影响,如果大就要对入口到圆柱之间的距离加大,或者适当调整一下网格。还有你看看你的courant number,我感觉这个不要太大,同样网格下,不同的courant number结果也会有不同。
-
piD 在昨天跑完10个周期之后,发现和piD/2,没有太大差别,1.3附近,稍小一点。
Courant number在整个计算中保持在0.2以下。
对于压差问题:
- 出口边界条件是 fixValue, 数值是0;
2.入口的压力边界条件设置是 zeroGradient;
流场还在重构中,传回我电脑还需要一点时间。等结果出来之后,我会回复你的。
谢谢你的回复。
- 出口边界条件是 fixValue, 数值是0;
-
@random_ran 不行,你就在试一试换一种网格画法,调整一下你的网格。
-
对于这个网格,个人感觉是能接受的。主要原因是用同样的网格用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的子网格上吗?
-
不好意思,是我犯错了。 fvScheme中用的是二阶精度的backward。
-
的确,在从o-shape形状区域到wake区域的转换有一点瑕疵。转换比太大?我会重画一下网格再试试。
-
首先,祝大家新年快乐! 对一个月前的问题继续进行补充。
用了 dynamicKEq 的子网格,得到的阻力和升力系数如下:
计算还在继续,虽然现在只算了几个周期。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 历史如下:
。他们的TKE模型计算得到的Cd的r.m.s 值很小,和我的试算反差很大。查看fvSolution文件后,我怀疑是残差标准的设置可能会影响数值计算结果,所以会想办法提高标准后看结果怎么样。
另外一个比较头疼的问题是尾流的速度曲线:文献中出现的这根曲线,横轴是计算域中轴线(原点在圆柱中心);纵轴是进口速度方向上的速度大小,通常是非常平滑的,如下图:
而我的计算结果,请注意是黑色虚线那条:
虽然大体上类似,但是波动却非常明显,不知道这是什么原因造成的。我的采样点是从Cl的历史中,发现振荡稳定后,取了10个周期。 而上篇文献则取了150多个周期,这个问题的来源吗?或许大家有更好的见解?
-
-
@random_ran 你这次只改了流场吗?你之前矩形的流畅有没有把过度的地方改一改,看看结果有影响吗?相同雷诺数,不同的流畅性质应该不会产生影响。
-
@sorrby 改了很多地方,比如fvm 的scheme, 时间步长,收敛标准。
对于之前那套网格,实在是没有精力了。一个二维的网格就16W的cell了,已经做得很密了。 反倒是圆形的这个计算域,只用到4W的cell就把很多参数算得很准了。要是有时间,再试试那个矩形的计算网格。
另外,我的意思是这样的,对于相同的雷诺数,不同的流体性质,会不会对一些测量参数产生一定范围内的影响?比如10%以内。
-
赞 @random_ran 做科研的态度:laughing:
-
-
Hi yhdthu:
我的这个case没有重点关注涡量的绝对大小。主要考察的是平均Cd, Strouhal 数,recirculation的长度。 voriticiy 只是用三维等值面做了一下wake的结构。 不知道你的case 重点考察的对象是什么? 如果你能提供的一下你的涡量考察的座标,我可以把我的数据提供给你。
对于蓝色的线的Cd 我是把每个一个时间步长的Cd都输出来。这个Cd是平均的整个Cylinder的表面。 然后把每个时间点的Cd 对时间用matlab plot一下就是那张图了。
-
@random_ran
好的,谢谢你
第一个应该是用的DES模型计算的吧,我看帖子最上面的Cd随时间变化幅度很小,但是后面哪个用了动力模式做的Cd图波动很大,虽然均值基本上是和实验对得上的我的目的是想做空化模拟,但事前要先将湍流结构捕捉的较好才可以,目前二维试算用的是动力模式,Cd波动也很大,我不是每个时间步输出的结果,平均起来效果不是很好,可能是样本点太少了
可否将Cp曲线也po出来呢?还有其时间平均的定义,希望能做个参考,非常感谢
BTW,如果可以的话,加个微信好友,可以随时交流~358253794
-
最开始确实是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:
-
@random_ran 好的,方法基本一致,我自己再整理一下,我觉得我的问题应该是后处理的问题,有问题再交流,祝好
-
近壁的地方是最难处理的,因为在这些地方没有大涡。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+计算第一层网格高度的公式。
可以看看这个帖子。
-
方形计算区,o型网格覆盖3D~4D范围,再加上尾迹区加密,结果应该是一样的。主要原因,猜测是非正交网格对涡解析不够,造成Cd或Cl差异。至于涡对Cd和Cl影响,可以参考吴介之老师的涡动力学。
-
当时的计算结果高了大概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网格的话,对计算结果影响大吗?
-
如果单纯从均值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_ran 在 Cd计算不准的问题-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中的修正。
-
主要参考的是
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是这么说的。可能更加有物理意义吧。 -
模型与模型之间可能是不一样的,参考NASA TMR 关于SA模型的实现,https://turbmodels.larc.nasa.gov/spalart.html#sa,SA模型有10种左右变种。
OF实现的是SA-fv3.
-
我觉得你第一个图做得已经很棒了!比对的是两组实验数据,你可以再加上数值模拟方面的数据就更好了。
如果你想再提高精度,更精密的网格肯定逃不了。另外湍流的模型的选择对结果的影响很大。还有就对流项的离散选择。关键是看你模拟的目的是什么。
Ps: 文献纵座标是废的是什么情况?能详细说说吗?
Ps2: 网格数量是多少?求解核心数?计算时间?边界条件如何?不知能否分享一下你的经验? -
@random_ran 你好,我的意思是纵坐标其实没有实际意义,可以想象,在y/D>=3时,Ux=Uinf,所以文献给出的其实是示意图,是为了方便说事。
另外,我用的就是方形计算域,大概是220W,壁面附近特殊处理了一下,分成三个增长率生成网格,过渡区网格注意了一下。基于升力的St大概是0.202,Cd约为1.02,下次在电脑旁我会把其他参数加上。
对流项离散LUST,就是那个75%CD,25%UD
DES类湍流模型,DDES_SA -
wallShearStress可以在计算完成后处理,O.F. (v4.1)是有这个功能的,但是个人觉得这个没法找到分离点的座标,只是给出了应力的最大和最小值。
simpleFoam -postProcess -func wallShearStress
在分离的地方,圆柱表面的剪切应力应该很接近零。不过一个点的剪切应力应该有三个分量,是考虑某个分量上的最小值还是合成的最小值?我也在想这个问题。
不过wallShearStress是和壁面速度沿壁面法向的梯度,所以用这个梯度也可以来求分离点的位置。结合ParaView是很容易实现的。
我建议还是用ParaView 后处理。
Filters-> Alphabetic->Gradient of Unstructured Dataset 对速度场进行处理。
然后用 Plot on intersection curve 输出圆柱表面的梯度分布,再输出所有时间步长的数据,就可以求出时间上的均值。如果在多切几个不同的界面,就可是实现空间上的均值。
希望能帮助到你。
发布 123 的 24