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

wwzhao

@yhdthu

被你搞晕了。。。说一下我是怎么分离点的吧。

修改wallGradU的代码，让它读取UMean，写出wallGradUMean。
得到wallGradUMean后，用paraview对圆柱表面的patch做slice，提取圆周上的点的wallGradUMean的值，保存到一个csv文件里。
用python写个脚本，读取并处理csv里面的数据，将wallGradUMean投影到圆周表面的切向方向。投影后的值出现正负变化的时候即为分离点。

yhdthu

@wwzhao 好的，我明白了，切向量是按着圆的切线方程得来的，这样上半圆与下半圆的切线方向是不同的吧？为了做出精确的曲线，在展向做切片的数量比较多如何处理？

wwzhao

@yhdthu

圆周上下是对称的，上下切一半即可。
对于展向切片较多的情况，建议你写个paraview的python脚本实现数据批量提取。

yhdthu

@wwzhao 你好我用你的方法试了，但是发现一个严重的问题，用paraview导出后的csv格式出现了数据丢失，比如我模型的x范围是(-1,1)，导出来之后区间缩小了，变成(-0.85,0.88)，请问这怎么处理呢？

wwzhao

@yhdthu 我没有遇到这个情况，你是怎么导出的？先做clip再做slice？

yhdthu

@wwzhao 我都试过，结果一样的，不知道怎么回事

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

CFD中文网

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