yPlus在openfoam代码里面的实现

xpqiu

u.boundaryField()[patchi].snGrad()

u.boundaryField()[patchi].snGrad()

加上 mag 试试，也就是改成

mag(u.boundaryField()[patchi].snGrad())

fangyuanaza

@xpqiu 谢谢老师~按您这样修改以后，确实没有sqrt的错误了，但是还有错误没有解决：

#0  Foam::error::printStack(Foam::Ostream&) at ??:?
#1  Foam::sigSegv::sigHandler(int) at ??:?
#2  ? in "/lib/x86_64-linux-gnu/libc.so.6"
#3  Foam::RASModels::kOmegaSSTyuPlus<Foam::EddyDiffusivity<Foam::ThermalDiffusivity<Foam::CompressibleTurbulenceModel<Foam::fluidThermo> > > >::correct() at ??:?
#4  ? in "/home/dyfluid/OpenFOAM/OpenFOAM-7/platforms/linux64GccDPInt32Opt/bin/rhoSimpleFoam"
#5  __libc_start_main in "/lib/x86_64-linux-gnu/libc.so.6"
#6  ? in "/home/dyfluid/OpenFOAM/OpenFOAM-7/platforms/linux64GccDPInt32Opt/bin/rhoSimpleFoam"
Segmentation fault (core dumped)

#3应该是最关键的，但是学生不太理解的是编译的时候仅仅是在基本的kOmegaSST湍流模型修改的BasicTurbulencemodel，这里为什么会出现fluidThermo名称空间下的错误?

xpqiu

@fangyuanaza

    forAll(uTau,patchi)

这里有问题，按照这样遍历，你其实遍历的是每一个体网格，假设你有1000个体网格，那么 patchi 的范围是0-999，但是你接下来的代码是把 patchi 当成边界的编号来用，肯定不行的。至于为什么出现 fluidThermo下面的报错，那个是因为湍流模型有着一个复杂的继承关系，可压缩的计算所使用的湍流模型就会用到 fluidThermo 类相关的东西。

fangyuanaza

@xpqiu 理解了，因此我应该取uTau的边界网格进行遍历

    forAll(uTau.boundaryFieldRef(),patchi)

或者

    forAll(uTau.boundaryField(),patchi)

可以运行成功！但是这里对这种编程思路不太理解：
由李老师提供的公式来看，这里
uTau应该是一个体标量场，这个值应该在近壁面是较大的，因为速度梯度较大，也就是边界层内部值比较大。那么为什么仅给网格边界赋值呢？而且，壁面上nut为零。如果只对边界进行赋值的话，uTau在壁面也是零。我输出uTau对比过了，uTau模拟出来结果在远场上边界才有值，其余地方为零，这是因为壁面nut为零，而场内部初始话赋予的零为初始值。

所以学生认为应该直接对整个网格场赋值。不知道理解对不对，请老师指点一下。

fangyuanaza

@xpqiu 请教老师，为什么不写成这样呢：

    forAll(uTau_, celli){
        uTau_[celli] = sqrt(nut[celli] * mag(tgradU()[celli].yx()));
    }
    uadd_ = u/(uTau_ + uTausmall);
    yadd_ = y_*uTau_/((this->mu()/this->rho_));

此处uTau定义为一个标量场。
是学生对这个公式理解不对吗？

xpqiu

@fangyuanaza
$u_{\tau}$ 应该是定义在壁面上的一个量，因为这个量跟壁面剪应力紧密关联，内部网格是没有壁面剪应力的。
首先，壁面剪应力的定义为
$$
\tau_w = \nu\frac{\partial u}{\partial y}
$$
$\nu$ 是分子粘度。
使用有限体积方法进行计算时，所有的计算都是离散的，因此
$$
\frac{\partial u}{\partial y} \approx \frac{u_{internal}-u_w}{\Delta y}
$$
$u_{internal}$ 表示第一层网格的速度，$u_w$ 是壁面速度，一般是0。
如果壁面附近网格画得很密（yPlus < 1），那么 $\frac{\partial u}{\partial y}$ 与 $\frac{u_{internal} - u_w}{\Delta y}$ 足够接近，因此这种情况下，壁面剪应力可以用
$$
\tau_w = \nu \frac{u_{internal}-u_w}{\Delta y}
$$

但是如果壁面网格画得粗，那么 $\frac{\partial u}{\partial y}$ 与 $\frac{u_{internal} - u_w}{\Delta y}$ 之间的差异就会很大。 OpenFOAM 的做法是利用 $\nu_t$ 来修正分子粘度，从而使得
$$
\tau_w \approx (\nu+\nu_t) \frac{u_{internal}-u_w}{\Delta y}
$$

回过头来看，其实不管是近壁网格画得粗还是细，都可以用上一个公式来计算 $\tau_w$ ，只不过壁面网格很细（ yPlus<1)时，可以认为 $\nu_t = 0$（但是 $\nu$ 不会等于0，所以壁面剪应力不可能为0，$u_{\tau}$ 也不可能为0）。而对于比较粗的网格，$\nu_t$ 比 $\nu$ 大得多，所以也可以认为 $(\nu_t + \nu) \approx \nu_t$，这应该也就是李老师写的公式的来源。

fangyuanaza

@xpqiu 感谢老师这么详细的指导！我这个算例是二维平板流动的，第一层y+<1，在粘性底层，所以nut为零，应该用nu。所以其实我根本不需要编程计算uTau，因为它可以由壁面监测的壁面摩擦力除以密度然后开根号得到。因为对于u+， y+曲线图来说，uTau在壁面固定地方是固定的。我也对比了这两个值，在边界上编程用nu代替nut，壁面摩擦速度数值是2.5657. 而用监测的壁面摩擦力及密度计算的壁面摩擦速度值是2.5734. 非常接近。再次感谢老师的指导，也感觉到看书的时候认为理解的概念，实际操作发现其实根本没懂，还是要多实践啊。谢谢老师~

vbcwl

@李东岳东岳老师我在做和这个前辈类似的东西，请问您有他的论文或者算例文件可以给我学习下么？

vbcwl

@xpqiu 您好，请问下我想利用utau来对距壁面距离和流体速度做无量纲化的，utau=sqrt(nuu0/y0) u0是壁面上方第一个网格的速度，y0是writeCellCentre命令写出的第一个网格的中心坐标，然后由于yPlus小于1，直接用mu来求解壁面剪应力，yPlus=utauy/nu,这样子算出来的yPlus与postProcess写出的yPlus有较大差距是什么原因？

xpqiu

@vbcwl
请问你计算的时候，nut 的壁面边界条件设置的是什么？
postProcess 算出来的 yPlus 是调用 nut 的边界条件里面定义的 yPlus 函数算出来的。

vbcwl

@xpqiu 感谢前辈的回复，我用的laminar做qDNS，没有用到nut这些，只用到了p p_rgh U T。请问您了解过没有nut文件的情况么？

xpqiu

@vbcwl
哦，如果是这样的话，那只能通过你描述的
$$
y^+ = \frac{y u_\tau}{\nu}
$$
的方式来算 $y^+$ 了。postProcess 这种计算方式调用的是湍流模型的东西，适合用湍流模型的场景。

shiyu

@Calf-Z-DNS 你好，看到你提的这个问题，请问你解决了吗？
我也想计算整个流动区域近壁面好几层网格的yPlus，而不是仅仅贴壁面的一层网格的yPlus，暂时我看到的wallShearStress和yPlus的计算都是贴壁面一层的值；
或者像东岳老师的示例代码那样直接给定了不变的uTau，而其实2维流动中，沿着流线方向uTau是变化的，所以我觉得比较复杂的是，怎么用某一壁面网格上的wallShearStress（或uTau）确定对应法线方向上的多层网格yPlus，这个对应关系不知道如何寻找。
希望我表达清楚了，非常感谢！

CFD中文网

yPlus在openfoam代码里面的实现