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

random_ran

再仔细分析了一下计算，代码，和计算结果，又产生了疑惑。

疑惑的对象是：二阶张量点乘(内积,dot product)一个矢量，其结果应该还是一个矢量。这个点乘与两个矢量之间的点乘不太一样。后者结果是一个数。其物理意义是一个矢量在另外一个矢量方向上的投影长度乘以被投影到的这个矢量的模长。

那么这里的二阶张量点乘一个矢量得到的这个新矢量代表着什么呢？

这个新的矢量是否代表这着 du/dn, dv/dn 以及dw/dn?

那这个新矢量的物理意义就是某个点的速度沿法向量方向上的梯度？并不是我们要求的速度沿圆周面切线的分量沿法向量的梯度？

越绕越糊涂。

wwzhao

@random_ran 之前我的表述有误。这里不涉及到二阶张量，snGrad返回的是一个矢量，代表壁面处速度沿法向的变化率。由于壁面上的速度为0，所以snGrad返回的实际上就是第一层网格的速度除以第一层网格格心到壁面的距离。snGrad是平行壁面的，在二维情况下（+x为右，+y为上，流速指向+x），取圆周上半部分，当snGrad的方向由顺时针变为逆时针时（snGrad在顺时针切向方向的投影符号发生变化），即为流动分离点。

wwzhao

@yhdthu 你看的是"Points:0"，"Points:1"和"Points:2"吗？这个不是按大小排序的，仔细检查一下。

yhdthu

@wwzhao 你好，我是直接把点转化成了角度值，发现只有60度到120度

yhdthu

@wwzhao 你好，我有个新思路，是直接在of代码中改的，具体如下：

想办法求取每点的切向量，其方向在上下半圆都是流动的反向，之后与以UMean计算的snGrad做内积

实现这个切向量的方法很简单，就是用一个vector(0,0,1)与面法向量做叉积，再乘以面法向量在y轴的投影正负号函数

这是最后的效果图：

wwzhao

@yhdthu

这个是我的结果，纵坐标已经处理过，分离角85.8度。

random_ran

再顺思路：

将速度(矢量)沿切线(矢量)分解，求出切线速度(矢量)
对1所得切线速度（矢量）求梯度，得到含九个分量的二阶张量
用这个张量点乘法向量得到一个新的矢量，也就是切线速度沿法线上的梯度
将3中所得矢量点乘切线向量，最后得到一个标量
沿圆周画出4得到的标量，找到等于零的点对应的座标

这是目前基于O.F.直接计算结果，我能想到的计算分离角的方法。ParaView似乎对叉乘不太支持？这样第一步似乎就走不通了。还需要再研究。

看来不动O.F.代码真是麻烦了...

yhdthu

@random_ran 我改的代码如下，写的可能有点啰嗦，如果有更好的表达形式望告知

        forAll(wallGradUMean.boundaryField(), patchi)
        {
            const fvPatch& currPatch = patches[patchi];
            if (isA<wallFvPatch>(currPatch))
            {
                	
                wallGradUMean.boundaryField()[patchi] =
                    (
                    	-UMean.boundaryField()[patchi].snGrad() &
                    	(
                    		(
            			mesh.Sf().boundaryField()[patchi]
            			/mesh.magSf().boundaryField()[patchi]
                    		) ^ vector(0,0,1)
                    	)
                    ) * sign(mesh.Sf().boundaryField()[patchi] & vector(0,1,0));
            }
        }

random_ran

一种方法：

求出平均速度之后，用ParaView的 GradientOfUnstructuredDataSet 直接计算vorticity在Z方向上的分量. 不用动代码，只需要ParaView 5.2.0 版本以上。

yhdthu

@random_ran 这个我也试过，还要再次处理，我觉得有点麻烦，所以改了代码，我比较懒😂

bingningmeng45

@random_ran 这个圆形正交网格正交性确实很好，我也是做一个圆柱绕流DNS的算例验证Re=3900，刚开始采用o型网格感觉网格质量不是很好，计算效率不是很高。想换成此篇文献中的圆形正交网格，但不知道这样的网格边界条件是怎么设定的，比如说出口，入口的边界是半圆吗？希望你能解答一下疑惑

李东岳

0_1499871510271_2017-07-12 22-57-27屏幕截图.png

最近在准备课程的算例，有一个就是圆柱绕流，我找了找文献，目前圆柱绕流期刊里面大部分都是这种圆形网格。

不知谁那有文献做过矩形网格和圆形网格的对比研究？

@random_ran
你那个更换计算网格之后好像结果明显提高很多？

bingningmeng45

@李东岳你这个边界条件怎么给的

李东岳

@bingningmeng45
固定速度进口，零法向梯度出口，
零法向梯度压力进口，固定值压力出口。
Re=150
我参考的是一个SCI的文章，不过那个文章里面并没有解释为什么要用圆形计算域。

如果单纯只是为了保证网格非正交质量，在圆形计算域外拓展为矩形不知有何影响？

http://www.sciencedirect.com/science/article/pii/S0307904X08000243

bingningmeng45

@李东岳 inlet区域是圆形网格哪个部分？左边半圆么？

李东岳

是的，左边半圆进口，右边半圆出口。

random_ran

@bingningmeng45
一开始的时候，我是用ICEM画的网格，矩形计算域。费劲周折对于3900的Re，Cd依然给出了偏高的预测。

圆形网格我是用ANSYS自带的网格生成，几乎是自动生成的。一方面质量高，另一方面操作简单，所以就没在纠结在矩形计算域为什么预测出偏高的Cd。 CFD中很重要的一步：网格无关性验证。就是说无论我用什么样形态的计算域，只要计算域足够大，随着对网格的加密，最终的数值解应该会去逼近一个值。很明显，用圆形计算域更容易实现这个过程。

我注意到你是用DNS来算。我不知道你的计算效率不高具体不高到什么程度。你用的什么版本的O.F.，你的计算设备怎么样？至于边界条件的设置，你可以参考摩托车算例。

cd $FOAM_TUTORIALST/incompressible/pisoFoam/les/

你可以试试

fixedValue

或者

turbulentInlet

我也用的是 @李东岳给出的解释。

@李东岳文献中这种圆形计算域越来越多。而且从圆形计算域衍生出来的计算域形态，比如在尾流曲渐进加密，可能能更好地捕捉尾流细节。

目前，我没有看到系统的比较计算域之间的差异，对这个经典问题的研究。还希望有了解的朋友指出。

这个网站上有人给出了这个经典算例，用网站的自己服务器两天左右就计算完成。

供你参考。

李东岳

@random_ran

目前，我没有看到系统的比较计算域之间的差异，对这个经典问题的研究。还希望有了解的朋友指出。

支持，可能专门做圆柱绕流的SCI里面或许有。

用网站的自己服务器两天左右就计算完成

你购买了他们的服务？2天大约用多少钱？

bingningmeng45

@random_ran foam-extend/3.1,广州天河2，用了864核算的，10s左右才能更新一次时间步，deltaT=0.001，效率不是很高，我画的是矩形区域O型网格，网格数量5千多万。计算300s，需要50多天时间，资源消耗太多，看别人的文献DNS圆柱绕流Re=3900资源消耗并不是很大,不知问题出在哪里；我现在准备换成圆形区域网格计算

random_ran

@李东岳
我是免费用户，那个网站很赞的地方是很多用户把模型都共享了。

@bingningmeng45
Hi，多谢分享细节，这个计算套装很酷。感觉你的网格可能画得太细了。

稍稍满足满足y+的要求就可以。

我也更倾向用圆形计算域。

另外可以试试 OpenFOAM-dev。大规模的计算，过多的processor文件会让你头痛的。

另外分解计算域的方式，不同scheme的选择，不同求解器的选择，不同的硬件，不同的节点选择,... 太多东西可以影响计算时间。

目前，我看到有人专门研究这些的是这篇文献：
S. Keough,2014,"Optimising the Parallelisation of OpenFOAM Simulations"。

但是计算量也不是很大，供你参考。

CFD中文网

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