nut壁面函数如何影响湍流模拟

chszkc

本人最近在学习OpenFOAM中的壁面函数。因为自己涉及到更多的是大涡模拟的计算，所以看的基本上是关于湍流运动粘度nut相关的内容。大概把理论部分过了一遍然后好奇OpenFOAM具体是如何执行的，就在几乎零C++认识的基础上看了下OpenFOAM的代码，发现有特别多不懂的地方。

一开始自己的疑问是，在大涡模拟中并没有与nut相关的方程，那么经过壁面函数计算得到的nut是如何进入到方程时间推进的计算中的？
看了下自己用的最多的 one-equation sgs tke 模型 kEqn的代码correct()

void kEqn<BasicTurbulenceModel>::correct()
{
    if (!this->turbulence_)
    {
        return;
    }

    // Local references
    const alphaField& alpha = this->alpha_;
    const rhoField& rho = this->rho_;
    const surfaceScalarField& alphaRhoPhi = this->alphaRhoPhi_;
    const volVectorField& U = this->U_;
    volScalarField& nut = this->nut_;         //nut在这里定义
    fv::options& fvOptions(fv::options::New(this->mesh_));

    LESeddyViscosity<BasicTurbulenceModel>::correct();

    volScalarField divU(fvc::div(fvc::absolute(this->phi(), U)));

    tmp<volTensorField> tgradU(fvc::grad(U));
    volScalarField G(this->GName(), nut*(tgradU() && dev(twoSymm(tgradU())))); //在这里用到了nut，并作为计算源项G的一个输入参数
    tgradU.clear();

    tmp<fvScalarMatrix> kEqn
    (
        fvm::ddt(alpha, rho, k_)
      + fvm::div(alphaRhoPhi, k_)
      - fvm::laplacian(alpha*rho*DkEff(), k_)
     ==
        alpha*rho*G  //生成项G在这里被使用
      - fvm::SuSp((2.0/3.0)*alpha*rho*divU, k_)
      - fvm::Sp(this->Ce_*alpha*rho*sqrt(k_)/this->delta(), k_)
      + kSource()
      + fvOptions(alpha, rho, k_)
    );

    kEqn.ref().relax();
    fvOptions.constrain(kEqn.ref());
    solve(kEqn);
    fvOptions.correct(k_);
    bound(k_, this->kMin_);

    correctNut();  // 求解sgs k 方程结束后，nut被更新，函数结束
}



void kEqn<BasicTurbulenceModel>::correctNut()
{
    this->nut_ = Ck_*sqrt(k_)*this->delta(); // 新一步得到的sgs k用来计算cell center上的nut
    this->nut_.correctBoundaryConditions();  // 基于wall function更新边界上nut的数值
    fv::options::New(this->mesh_).correct(this->nut_);

    BasicTurbulenceModel::correctNut();
}

由上面的代码看来，nut在边界处的数值似乎并没有进入到矩阵求解的系数中（个人的理解，如有不对烦请指正），但从下面这段代码可以看到，calcNut()是在更新壁面上的nut数值，而并没有修改第一层网格中心上nut的数值

tmp<scalarField> nutkWallFunctionFvPatchScalarField::calcNut() const
41 {
42 const label patchi = patch().index();
43
44 const turbulenceModel& turbModel = db().lookupObject<turbulenceModel>
45 (
46 IOobject::groupName
47 (
48 turbulenceModel::propertiesName,
49 internalField().group()
50 )
51 );
52
53 const scalarField& y = turbModel.y()[patchi];
54 const tmp<volScalarField> tk = turbModel.k();
55 const volScalarField& k = tk();
56 const tmp<scalarField> tnuw = turbModel.nu(patchi);
57 const scalarField& nuw = tnuw();
58
59 const scalar Cmu25 = pow025(Cmu_);
60
61 tmp<scalarField> tnutw(new scalarField(patch().size(), 0.0));
62 scalarField& nutw = tnutw.ref();
63
64 forAll(nutw, facei)
65 {
66 label faceCelli = patch().faceCells()[facei];
67
68 scalar yPlus = Cmu25*y[facei]*sqrt(k[faceCelli])/nuw[facei];
69
70 if (yPlus > yPlusLam_)
71 {
72 nutw[facei] = nuw[facei]*(yPlus*kappa_/log(E_*yPlus) - 1.0);
73 }
74 }
75
76 return tnutw;
77 }

我个人的理解是，可能在经过壁面函数更新后，boundary上的nut经过某段代码赋值给第一层网格中心点上的nut，使得经过壁面函数修正的nut得以进入迭代计算中。可惜的是本人暂时还没找到执行上述猜想的代码，那表明OpenFOAM实际的执行方式与我理解的不同？

知乎上一位前辈似乎也提到了相同的疑问
https://zhuanlan.zhihu.com/p/32520364

这么说我的猜测的对的咯？那为什么我copy其中一个值，在internal field里搜不到？难道数据结构跟我想的不一样？那internal field里对应点上的nut值又是什么？当计算到最靠近壁面的那个cell时，到底应该用internal field里的nut还是boundary上的nut？

我目前的猜测，对于最靠近壁面那个点，其实有两个值，一个值是根据正常turbulence model算出来的，储存在该变量的internal field里；另一个是根据wall function算出来的，储存在该变量的boundary field里。然后真正计算矩阵的时候，会检查一下有没有用wall function，假如用了，就忽略internal field里的值，把wall function算出来的值覆盖上去。这里没有直接填充是为了数据处理安全，反正数据多了可以不用，少了就麻烦了。

但好像问题还是没有被解决，毕竟没有贴出来相应模块的代码......想请各位对这方面有深刻理解OFer指教指教，非常感谢！

李东岳

这部分内容还是值得讨论的。明天我放在《笔记》里面，更新后告诉你

李东岳

@chszkc 在 nut壁面函数如何影响湍流模拟中说：

在大涡模拟中并没有与nut相关的方程，那么经过壁面函数计算得到的nut是如何进入到方程时间推进的计算中的？

$\nu_t$在动量方程中起作用。在这一项里面$\nabla\cdot((\nu+\nu_t)\bfU\bfU)$，不在湍流模型里面。湍流模型只是求解。

nut在边界处的数值似乎并没有进入到矩阵求解的系数中（个人的理解，如有不对烦请指正）

$\nu_t$不存在PDE方程，因此你说的是对的。

boundary上的nut经过某段代码赋值给第一层网格中心点上的nut，使得经过壁面函数修正的nut得以进入迭代计算中。

不是。
网格第一层的$\nu_t$是通过$\varepsilon,k$来求解的，跟边界没有关系。外链网站的内容太多，没细看，暂不回答咯

总体来讲，壁面函数的流程主要是手动计算壁面第一层的epsilon以及湍流产生项（在OpenFOAM里面一般称之为G），然后$\nu_t$通过代数公式求出即可。更详细的我更新在这里了 http://www.dyfluid.cn/theory.pdf http://dyfluid.com/docs/wallFunction.html

chszkc

@李东岳非常感谢！我也稍微想通一点了，壁面函数的作用只会出现在边界处的nut上，时间步第n步通过壁面函数的到的patch上的nut，会在第n+1步的动量方程粘性项求解中起作用所以第一层网格中心上的nut只要通过k求解得到就好了，不需要被壁面函数更新。希望这次没有理解错了