请问关于wallShearStress 和 grad(U)

sibo

大家好，

最近在看wallShearStress的代码，发现OpenFOAM的计算是：

   forAll(wallShearStress.boundaryField(), patchI)
   {
       wallShearStress.boundaryField()[patchI] =
       (
          -mesh.Sf().boundaryField()[patchI]
          /mesh.magSf().boundaryField()[patchI]
       ) & Reff.boundaryField()[patchI];
   }

这里的Reff发现是：

const volSymmTensorField Reff(model->devRhoReff());

devRhoReff() 返回的应该就是 stress tensor.
那么这样的话，openfoam计算的wallshearstress应该是traction而不是shear stress啊.
大家可以参考下这个链接：traction
请问是这样吗？

后来想用wallGrad(U)来计算shear stress，发现OF5没有这个命令了。只能用

solverName -postprocess -func grad(U)

这里我有两个问题，
第一个是solverName对结果有关系么？因为我是用的自己的solver，和rhoCentralFoam类似，没有用simple或者PISO算法，所以这里如果用自己solver的名字，会报错。用rhoCentralFoam可以运行。按我理解只要有U这个场，有wall function 就应该可以可以了。所以想问大家，这个name对结果有影响么？

第二个问题是，用grad(U)后生成是是一个tensor field，每一个cell会有9个分量。但是他们对应的是依次怎么排序的呢？
是dU_x/dx, dU_x/dy, dU_x/dz, dU_y/dx, dU_y/dy, dU_y/dz, dU_z/dx, dU_z/dy, dU_z/dz ? 还是dU_x/dx, dU_y/dx, dU_z/dx... 呢？

非常感谢大家的帮忙！

李东岳

第一个是solverName对结果有关系么？

有的处理器没有grad(U)这个功能（没有包含相应的库或者场数据不够），这只能在尝试之后发现。你用自己的求解器，需要把后处理这面包含，同时，需要创建相应的场。

第二个问题是，用grad(U)后生成是是一个tensor field，每一个cell会有9个分量。但是他们对应的是依次怎么排序的呢？是dU_x/dx, dU_x/dy, dU_x/dz, dU_y/dx, dU_y/dy, dU_y/dz, dU_z/dx, dU_z/dy, dU_z/dz ? 还是dU_x/dx, dU_y/dx, dU_z/dx… 呢？

这个问题我在好多期刊都看到不同的答案。有的是转置了，有的没转。我还专门发了个朋友圈。我觉得是这么排序的（CFD这面），公式6。你觉得呢？

或许直接看看grad()函数的代码更好

0_1537409786791_微信图片_20180920101600.jpg

sibo

@东岳感谢李老师回复！
其实我是打算把grad(U)的plot和wallShearStress的plot对比一下看那个分量可以使他们有相似的shape哈哈。因为按理说同一个boundary上面的grad(U)和shear stress应该是差了一个mu(viscosity)的倍数。
0_1537411816267_072385a8-aac1-4cd1-9252-2f423172e99c-image.png https://www.cfd-online.com/W/images/math/0/c/b/0cbe8a602784d6e0c0c532d495dddff5.png

然而看起来所有的分量都对不上。。。难道是对这个物理量理解错误么还是说openfoam计算shear stress或者grad(U)的方法比较特别？？不过的确直接用wallShearStress command 导出的是traction 而不是真正的shear stress。就像我上个post里说的那样，不知道东岳师兄觉得是不是？

谢谢！

random_ran

1 应力

我看了你给出的参考链接， Traction 的定义是应力在边界垂直方向上的分量 (Wolfgang)。也就是说 Traction 是一个矢量 (3x1), 而应力应该是一个2阶张量（3x3).

我也发现 devRhoReff 的返回值是 effective stress tensor.

但是为什么楼主认为 OF 计算的就是 traction 了呢？能再解释一下吗？

2 后处理的问题

在你求解器的主函数中，添加：

#include "postProcess.H"

就可以在你开发的求解器中用后处理了。

sibo

@random_ran
感谢您的回复！

根据那个连接，traction是shear stress的在边界y方向的分量。这里如果我们定义单位法向量用n表示，stress tensor用T表示。那么，根据OpenFOAM的代码，
```
n = -mesh.Sf().boundaryField()[patchI] / mesh.magSf().boundaryField()[patchI];

T = Reff.boundaryField()[patchI];
```
那么OpenFOAM这里做的就是 n * T 啊，也就是traction的定义。所以按照那个链接最终给出的公式，

正确的计算应该是：
wallShearStress = nT - n(n*T) *n
加好了！非常感谢！

random_ran

@sibo

这个链接是 OF 给出的墙剪力计算的描述。

这里 OF 明确给出墙面剪应力是应力张量 ( R ) 与垂直于墙面的向量( n ) 的内积。一开始我也认为（请原谅我有限的数学底子），要算墙剪力，内积的对象应该是切于墙的向量。但是我后来发现，这个 Stress ( R ) 有点不一样。按照 OF 给出 devReff 的定义： devReff 是 Deviatoric part of the effective Reynolds stress。而这个 Deviatoric part 则是这个问题的关键。 Wiki 对 deviatoric 的解释是任何一个应力张量都可以写成两个特殊的应力张量之和：

hydrostatic stress (平均应力，对角矩阵)
deviatoric stress (总应力矩阵中减去 hydrostatic stress)

而 OF 返回的就是第2个应力张量, 在 devRhoReff 的定义（不可压中devReff?）中（L84~L85）:

(-(this->alpha_*this->rho_*this->nuEff())) * dev(twoSymm(fvc::grad(this->U_)))

( 可以参考 twoSymm 和 dev )

有效动力黏度 ($\mu$) 乘以经过 deviatoric 变换的速度梯度张量，得到的就是“有效应力张量”。最后，用“有效应力张量”与垂直墙的单位向量作内积，就得到了墙剪力。

可问题是，为什么要绕这么大一圈，而不是用：

\begin{equation}
\tau = \mu \frac{\partial U_x}{\partial y}
\end{equation}

求出垂直于墙面的速度梯度然后乘以运动黏度这样直接的计算方式呢 (湍流边界层中的运动黏度不是常量)？还请大家指点。

wwzhao

受OpenFOAM中wall function实现方式的影响，贴近壁面处的du/dy并不是真实的速度梯度，因此用du/dy求出来的并不是真实的wall shear stress。

也就是说du/dy的算法只适用于laminar flow，而不适用于turbulent flow。