Cd计算不准的问题-LES-Re3900-三维圆柱绕流-pisoFoam
-
@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 输出圆柱表面的梯度分布,再输出所有时间步长的数据,就可以求出时间上的均值。如果在多切几个不同的界面,就可是实现空间上的均值。
希望能帮助到你。
-
希望讨论能把这些疑惑逐渐拨开。
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吧
-
-
- 我知道如何求法向量(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 这样看来,并不需要再乘以面法向向量了