Cd计算不准的问题-LES-Re3900-三维圆柱绕流-pisoFoam
-
希望讨论能把这些疑惑逐渐拨开。
wallShearStress和近壁wall-normal的梯度成比例关系,所以 du/dn=0 的点就是wallShearStress为零的点。公式du/dn=0 中u指得是平行于壁面的速度,而n是壁面指向流体的法向量。三维问题不好描述,先在二维平面讨论。
假如我们考察圆周上(此时我们在二维平面内)任意一点A。A点处的速度在x轴和y轴(直角座标,原点在圆心,3点钟方向是x轴正方向,12点方向是y轴正方向)必有两个分量。但我们关心的是A点处的切向速度的大小:
$$ \v_{\tau} = v \cdot \tau) $$ 。
有了切向速度之后,在A点对这个速度沿A点指向流体的法向量求偏导,就可以得出在A点处平行于壁面的速度沿该点处壁面法向量的偏导。 (感觉很拗口,不知道这样理解对不对)
然而实际用ParaView来计算的时候还是有疑惑。
- 我知道如何求法向量(normal vector),切向量怎么求?
- 假如我求出了切向量并得到了切向速度大小,是否就可以直接求出在A点处平行于壁面的速度沿该点处壁面法向量的偏导,而不是这个速度分别沿x,y的偏导?
最初的想法其实很简单,用 Filters-> Alphabetic->Gradient of Unstructured Dataset 对速度场进行处理。处理完后得到一个2阶张量,这个二阶张量没有直接给出在A点处平行于壁面的速度沿该点处壁面法向量的偏导。但是直觉告诉我,这个2阶张量经过一定的"变换"是可以得到在A点处平行于壁面的速度沿该点处壁面法向量的偏导。
越来越把自己绕进去了。寻求大家的帮助。
-
$\nabla \mathbf U \cdot \mathbf n$=grad(U)*normal吧
-
既然分离处mag(wallGradU)是0,那么每个时刻云图的标尺都应该有0才对
不是 mag(wallGradU) 为0,而是wallGradU在法向的分量为 0,你需要做一下坐标变换才行。
我用lambda2处理isosurface时,发现出来的图不对,不该有涡的地方却出来了很不合理的结构
看你后处理用的应该是Tecplot,这个软件在对OpenFOAM算例做Contour的时候,会在zeroGradient类型的边界上也出现isosurface,算是个bug吧。你可以试试Paraview。
-
- 我知道如何求法向量(normal vector),切向量怎么求?
法向量与切向量点乘等于0,这个很好求吧
- 假如我求出了切向量并得到了切向速度大小,是否就可以直接求出在A点处平行于壁面的速度沿该点处壁面法向量的偏导,而不是这个速度分别沿x,y的偏导?
不用这么麻烦,做个坐标变换就好了。将梯度从x-y坐标系变换到法向-切向坐标系下。
-
关于分离点的做法,重理一下思路。
流体在空间中运动,比较简单的情况下(比如沿A点流向B点)(不太严谨),速度会减小, A点的静压力会高于B点的静压(伯努利方程)。我们说 (p_B - p_A)> 0。 也就是正压力梯度。
对于边界层出现流体分离的问题,比如我们研究的这个圆柱绕流问题。 在边界层的一定区域,由于复杂的机理(暂时没想明白),出现了负压力梯度。这就意味着,在这个区域内的流体有向逆方向运动的趋势。从总体来看,流体是沿正方向流动。不过总会有一个瞬间一个时刻,这个逆潮流而上的部位,把这个部分的流体转向。在这个区域内,总有一个点的速度趋近与零。这个点就是分离点。
如何找到这个分离点呢?
现在把目光放在圆柱体表面,来研究沿圆柱表面法向量方向(指向流体)的速度分布情况。在没有分离出现的部位,我们沿法向量方向(指向流体)走一小步,速度会增大(从壁面表面出发的速度为零)。再考虑发生分离部分的流场,如果走同样的路径,你会发现速度同样会增大,但是速度的方向和之前正好是反向。 如果考虑临界情况,就是说,当我们走出这一小步的时候,发现速度竟然没有变化!这个极限,就是分离点!如果把我们沿法向量方向(指向流体)走一小步并考察速度的变化,换成数学语言,就是du/dn=0的点!
接下来我再说说如何求这个点。
数学上讲,就是求速度沿壁面法向量的梯度。求梯度在ParaView可以通过
Filters-> Alphabetic->Gradient of Unstructured Dataset 对速度场进行处理,得到的是一个含有9个分量的二阶张量。但是这个并不是我们要求的沿壁面法向量的梯度。不过,还记得点乘(内积)的物理意义吗?如果用这个二阶张量”点乘” 法向量这个向量,不就得到了梯度在法向量那个方向上的大小了吗?ParaView对速度场求梯度之后在每一个点得到了9个分量,从gradient 0 一直到 gradient 8。如果没有猜错,gradient 0 代表d(U_x)/dx, gradient 1 代表d(U_x)/dy, gradient 2 代表 d(U_x)/dz。 这样的话用这个含九个分量的二阶张量点乘法向量,就可以得到速度梯度在法向量方向上的大小。计算方法如下:
wallshearStress(成比例于速度沿法向量梯度,不是严格相等) = sqrt((Gradients_0*Normals_X + Gradients_1*Normals_Y + Gradients_2*Normals_Z )^2+(Gradients_3*Normals_X + Gradients_4*Normals_Y + Gradients_5*Normals_Z)^2+(Gradients_6*Normals_X + Gradients_7*Normals_Y + Gradients_8*Normals_Z)^2)
计算完成之后,得到的是一个标量。从这个图可以看出来过了极大值后,这个量开始骤降,第一次出现极小值点的地方,就应该是分离点了。
还请大家批评指正。
-
@random_ran 说得很好,我是通过修改wallGradU的代码实现的,直接给出了时间平均的值,不过在GradU为0的点后面的部分,需要人工手动加负号,体现速度方向相反,关键代码如下:
wallGradUMean.boundaryField()[patchi] = -UMean.boundaryField()[patchi].snGrad();
snGrand()就是面法向梯度方向了,显然这是个vector,mag一下即可
@wwzhao 这样看来,并不需要再乘以面法向向量了
-
再仔细分析了一下计算,代码,和计算结果,又产生了疑惑。
疑惑的对象是: 二阶张量点乘(内积,dot product)一个矢量,其结果应该还是一个矢量。这个点乘与两个矢量之间的点乘不太一样。后者结果是一个数。其物理意义是一个矢量在另外一个矢量方向上的投影长度乘以被投影到的这个矢量的模长。
那么这里的二阶张量点乘一个矢量得到的这个新矢量代表着什么呢?
这个新的矢量是否代表这着 du/dn, dv/dn 以及dw/dn?
那这个新矢量的物理意义就是某个点的速度沿法向量方向上的梯度?并不是我们要求的速度沿圆周面切线的分量沿法向量的梯度?
越绕越糊涂。
-
@random_ran 之前我的表述有误。这里不涉及到二阶张量,snGrad返回的是一个矢量,代表壁面处速度沿法向的变化率。由于壁面上的速度为0,所以snGrad返回的实际上就是第一层网格的速度除以第一层网格格心到壁面的距离。snGrad是平行壁面的,在二维情况下(+x为右,+y为上,流速指向+x),取圆周上半部分,当snGrad的方向由顺时针变为逆时针时(snGrad在顺时针切向方向的投影符号发生变化),即为流动分离点。
-
再顺思路:
-
将速度(矢量)沿切线(矢量)分解,求出切线速度(矢量)
-
对1所得切线速度(矢量)求梯度,得到含九个分量的二阶张量
-
用这个张量点乘法向量得到一个新的矢量,也就是切线速度沿法线上的梯度
-
将3中所得矢量点乘切线向量,最后得到一个标量
-
沿圆周画出4得到的标量,找到等于零的点对应的座标
这是目前基于O.F.直接计算结果,我能想到的计算分离角的方法。ParaView似乎对叉乘不太支持?这样第一步似乎就走不通了。还需要再研究。
看来不动O.F.代码真是麻烦了...
-
-
@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)); } }
-
一种方法:
求出平均速度之后,用ParaView的 GradientOfUnstructuredDataSet 直接计算vorticity在Z方向上的分量. 不用动代码,只需要ParaView 5.2.0 版本以上。
-
@random_ran 这个我也试过,还要再次处理,我觉得有点麻烦,所以改了代码,我比较懒😂
-
@random_ran 这个圆形正交网格正交性确实很好,我也是做一个圆柱绕流DNS的算例验证Re=3900,刚开始采用o型网格感觉网格质量不是很好,计算效率不是很高。想换成此篇文献中的圆形正交网格,但不知道这样的网格边界条件是怎么设定的,比如说出口,入口的边界是半圆吗?希望你能解答一下疑惑
-
@bingningmeng45
固定速度进口,零法向梯度出口,
零法向梯度压力进口,固定值压力出口。
Re=150
我参考的是一个SCI的文章,不过那个文章里面并没有解释为什么要用圆形计算域。如果单纯只是为了保证网格非正交质量,在圆形计算域外拓展为矩形不知有何影响?
http://www.sciencedirect.com/science/article/pii/S0307904X08000243
-
@bingningmeng45
一开始的时候,我是用ICEM画的网格,矩形计算域。费劲周折对于3900的Re,Cd依然给出了偏高的预测。圆形网格我是用ANSYS自带的网格生成,几乎是自动生成的。一方面质量高,另一方面操作简单,所以就没在纠结在矩形计算域为什么预测出偏高的Cd。 CFD中很重要的一步:网格无关性验证。就是说无论我用什么样形态的计算域,只要计算域足够大,随着对网格的加密,最终的数值解应该会去逼近一个值。 很明显,用圆形计算域更容易实现这个过程。
我注意到你是用DNS来算。我不知道你的计算效率不高具体不高到什么程度。你用的什么版本的O.F.,你的计算设备怎么样?至于边界条件的设置,你可以参考摩托车算例。
cd $FOAM_TUTORIALST/incompressible/pisoFoam/les/
你可以试试
fixedValue
或者
turbulentInlet
我也用的是 @李东岳 给出的解释。
@李东岳 文献中这种圆形计算域越来越多。而且从圆形计算域衍生出来的计算域形态,比如在尾流曲渐进加密,可能能更好地捕捉尾流细节。
目前,我没有看到系统的比较计算域之间的差异,对这个经典问题的研究。还希望有了解的朋友指出。
这个网站上有人给出了这个经典算例,用网站的自己服务器两天左右就计算完成。
供你参考。
-
目前,我没有看到系统的比较计算域之间的差异,对这个经典问题的研究。还希望有了解的朋友指出。
支持,可能专门做圆柱绕流的SCI里面或许有。
用网站的自己服务器两天左右就计算完成
你购买了他们的服务?2天大约用多少钱?
-
@random_ran foam-extend/3.1,广州天河2,用了864核算的,10s左右才能更新一次时间步,deltaT=0.001,效率不是很高,我画的是矩形区域O型网格,网格数量5千多万。计算300s,需要50多天时间,资源消耗太多,看别人的文献DNS圆柱绕流Re=3900资源消耗并不是很大,不知问题出在哪里;我现在准备换成圆形区域网格计算
-
@李东岳
我是免费用户,那个网站很赞的地方是很多用户把模型都共享了。@bingningmeng45
Hi,多谢分享细节,这个计算套装很酷。感觉你的网格可能画得太细了。稍稍满足满足y+的要求就可以。
我也更倾向用圆形计算域。
另外可以试试 OpenFOAM-dev。 大规模的计算,过多的processor文件会让你头痛的。
另外分解计算域的方式,不同scheme的选择,不同求解器的选择,不同的硬件,不同的节点选择,... 太多东西可以影响计算时间。
目前,我看到有人专门研究这些的是这篇文献:
S. Keough,2014,"Optimising the Parallelisation of OpenFOAM Simulations"。但是计算量也不是很大,供你参考。
-
@random_ran
Mean streamwise velocity along the centreline of the cylinder这幅图怎么进行z轴(圆柱轴向方向)平均?还有对比其他的量也需要z轴平均。需要在算例怎么设置 -
-
@random_ran 你这个做法我不是很懂,能详细说明一下吗,谢谢
-
就是把ParaView中的操作脚本化。然后对生成的脚本做一些“胶水”处理。我这里用的Perl,你可以用自己熟悉的编程语言。
获得*.py 脚本的内容这样操作,在ParaView GUI 里面找到:
Tools - Start Trace
进行你的图形操作,直到输出数据文件为止。
然后 Stop Trace,就生成了 *.py 的脚本。
接下来对这个脚本文件处理就好。这个脚本文件中绝大部分都是一样的,只要你在感兴趣的地方稍作修改就好,比如储存地址,切片位置。
用Perl的目的就是少了做了几次“复制”“粘贴”。
最后生成一个“真正”的脚本,执行
pvpython *.py
-
我现在正在运行一个圆柱绕流Re=3900的算例,在处理一维波数谱图(1D wavenumber spectra)中遇到问题,不太清楚下图横纵坐标如何处理得到。
纵坐标是u'u'么?
横坐标是波数k_z,如何得到的?大部分的文献说的很笼统,并不能很清楚知道他们的计算处理过程。
我现在的做法是沿圆柱轴线z方向,截取一条线上的流向(x轴)速度波动,横纵坐标不太会处理,一直纠结在这。 -
我从维基里面扣下来的一段话:
"从随着时间而变的函数萃取出的一组数据,经过傅里叶变换,会得到一个频率谱;而从随着位置而变的函数萃取出的一组数据,经过傅里叶变换,会得到一个波数谱。"
相信对于速度在时间历程上进行快速傅里叶变换(FFT),应该很熟悉。时间域内(时间, 某个方向上的速度大小) --> (频率,某个方向上的速度大小); 对应的 位置域内(位置,某个方向上的速度)--> (波数,某个方向上的速度大小)。 [位置域内是我瞎造的], 不过量纲上是讲的过去。 时间FFT后 横轴变换成 s-1; 位置(横座标)FFT后变成 m-1 也就是波数的单位。
我的感觉是作者取了4条直线。 每条直线上搜集的是沿流向方向上的速度的波动值 (u'u') [对应雷诺应力],然后以 (纵座标,u'u') 做FFT。
但是为什么这样做呢?
-
@random_ran 你好,请教nut和nuTilda边界怎么设置
-
@random_ran @yhdthu 我对于那条蓝色的cd曲线还存在问题,
这个Cd是平均的整个Cylinder的表面。
这个是如何平局得到的呢?
我的理解是计算输出每一步的结果,从某个时刻开始,开始对瞬时值做时间平均,是这样的吗? -
-
nut ($\nu_t$) 是湍流粘度;
-
nuTilda ($\overset{\sim}{\nu}$) 是 Spalart-Allmaras 模型中引入的一个量
对于 $\overset{\sim}{\nu}$:
- 墙上: 0
- 进口处: 0 最好
SA 模型认为 $\overset{\sim}{\nu}$ 小于 $\nu/10$的数值是可接受的。同样的情况适用于初始条件。
另外, 这个模型还强调了对 free-slip 边界条件的使用约束: Neumann (法向梯度为零)。 以上是参考:
Spalart, P.R., & Allmaras, S.R. (1994).
A one-equation turbulence model for aerodynamic flows.
La Recherche Aerospatiale, 1, 5-21.当然,最保险的方式是参考算例模型了。
@benqing 对的,Cd 是对每个时间步上的瞬时值进行时间上的平均来得到均值的。 注意,这个瞬时值已经是整个 Cylinder 表面的均值。
@东岳 投到一本书上,已经接受啦。
-
-
这是我的毕业论文,纪念这个帖子盖到 100 楼。
虽然已经毕业,但依旧是 CFD 的小白。每次打开主页的时候,大把大把的专业名词都是我从来没接触到的。但我依旧会去探索 CFD, 探索 OF。
一路走来,学习 OF 的路上,真的是很感谢那些曾经帮助过我的人。 特别是这里,‘ cfd-china.com’,很喜欢这个论坛,简洁的界面,却不失复杂的功能,当然更重要的是, 还有藏龙卧虎的你们 : )。
感谢你们。