yPlus在openfoam代码里面的实现

[李东岳]

最近算了几个验证壁面函书的LES，可能加van Driest会更准。有空把这个算立整理整理发上来

另外得知：

@东岳 在 yPlus在openfoam代码里面的实现 中说：

如何在非壁面处求yPlus的数学公式？

自我回答，跟常规方法一样。程序里：

#include "wallDist.H"

int main()
{
    volScalarField Dist(wallDist::New(mesh).y());
    dimensionedScalar uTau("uTau", dimVelocity, 1.003);
    volScalarField yPlus(Dist*uTau/turbulence->nu());
}

[Calf.Z-DNS]

@东岳 谢谢老师！这里面的uTau为何直接取1.003呢？不需要植入它的公式吗？

[Calf.Z-DNS]

如果要植入MK湍流模型的话，里面有一项
其中y+我按照了东岳老师提供的方法，那就是如下图：

编译成功，但按这个湍流模型跑的结果不对。

uTau为何设成1.003呢？如果要植入uTau的公式，那tauw这个如何得到呢？源代码里面相关的有tauw跟wallshearstress，但是不知道怎么引用。还望赐教。

[李东岳]

@Calf-Z-DNS 我那个是我自己写哒，计算出来之后是1.0003，具体我有点忘了下周回办公室给你看看

[Calf.Z-DNS]

@东岳 好的谢谢！

[Calf.Z-DNS]

@东岳 老师 ，请问uTau不应该用公式植入吗？在沿管方向的壁面上，uTau应该是变化的吧。如果用初始值计算出来的uTau定值植入，我试了下，湍流模型没有达到自己预想的效果。

[李东岳]

@Calf-Z-DNS 对，是用公式计算的，明天去办公室给你看看

[李东岳]

$u_\tau=\sqrt{\nu_t \frac{\rd u}{\rd y}}$

我在几年前看过有人做管道流动的温度DNS，非常酷炫，一直在关注你的结果，有好结果了发上来炫一波啊

[Calf.Z-DNS]

@东岳 谢谢老师。uTau的公式我知道，但不用把这个公式植入代码里面吗？如果需要的话，有点不太清楚怎么写入代码，还望赐教。
DNS圆管的算例目前还在算呢，所需要的资源太大了，算比较慢。

[李东岳]

@Calf-Z-DNS

forAll(xxx, patchi)
{
    utau.boundaryFieldRef()[patchi] = 
        sqrt(nut.boundaryField()[patchi]*u.boundaryField()[patchi].snGrad());
}

[李东岳]

我之前见过一个视频，DNS模拟这种平板之间的浮力驱动流，最终里面的流体呈现这种弯曲的流线，这叫什么现象？最近打算复现一下这个东西

[李东岳]

@东岳 https://arxiv.org/pdf/1904.11400.pdf

[Calf.Z-DNS]

我做的压力驱动流，也是两个平板间的channel flow，上下有温差，但不考虑浮力。瞬时温度场长这样：

[李东岳]

你这个有弯曲的形状出来么？

本来打算这次OKS准备下这个模拟时间不够了，下次OKS准备一下，https://arxiv.org/pdf/1904.11400.pdf 主要是这个工况里面的流动感觉很有意思

[Calf.Z-DNS]

没有具体看流线，但我觉得应该是弯曲的，不过不清楚是否是规律弯曲的。

[fangyuanaza]

@李东岳 请教下李老师~对于这里计算uTau，我用dimensionedScalar 是可以运算的，但是用这块代码计算uTau会报错，可能理解不太对。对于forAll中的×××，认为是uTau, 在这之前，我先定义了uTau是一个体标量场，初始值给的0，您代码中的u.boundary中的u认为是U.component(0)(), 这样代码能编译成功，但是运行几步会报错，请老师指教一下，万分感激(Ps 我在湍流模型里面实现这几个量，离壁面距离已经给出了，就是y_)：

    volScalarField u(U.component(0)());
    volScalarField uTau(0*sqrt(k_));
    //dimensionedScalar uTau("uTau", dimVelocity, 1.003);
    dimensionedScalar uTausmall(
       "0", dimensionSet(0, 1, -1, 0, 0, 0, 0), 1e-20
    );
    forAll(uTau,patchi)
    {
        uTau.boundaryFieldRef()[patchi] = sqrt(nut.boundaryField()[patchi]*u.boundaryField()[patchi].snGrad());
    }
    uadd_ = u/(uTau + uTausmall);
    yadd_ = y_*uTau/((this->mu()/this->rho_));

运行报错为：

#0  Foam::error::printStack(Foam::Ostream&) at ??:?
#1  Foam::sigFpe::sigHandler(int) at ??:?
#2  ? in "/lib/x86_64-linux-gnu/libc.so.6"
#3  Foam::sqrt(Foam::Field<double>&, Foam::UList<double> const&) at ??:?
#4  Foam::sqrt(Foam::tmp<Foam::Field<double> > const&) at ??:?
#5  Foam::RASModels::kOmegaSSTyuPlus<Foam::EddyDiffusivity<Foam::ThermalDiffusivity<Foam::CompressibleTurbulenceModel<Foam::fluidThermo> > > >::correct() at ??:?
#6  ? in "/home/dyfluid/OpenFOAM/OpenFOAM-7/platforms/linux64GccDPInt32Opt/bin/rhoSimpleFoam"
#7  __libc_start_main in "/lib/x86_64-linux-gnu/libc.so.6"
#8  ? in "/home/dyfluid/OpenFOAM/OpenFOAM-7/platforms/linux64GccDPInt32Opt/bin/rhoSimpleFoam"
Floating point exception (core dumped)

说是这个里面的sqrt传入变量不对？

[xpqiu]

@fangyuanaza 在 yPlus在openfoam代码里面的实现 中说：

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为零，而场内部初始话赋予的零为初始值。

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