请教各位老师~
****目标：***希望给雷诺应力增加高阶项来修改湍流模型，高阶项的形式是固定的（标量对称张量）
***想法：***因为雷诺应力本身就是对称张量，其中的标量是随时变化的，机器学习生成的。如果用include **.H的方式来修改湍流模型，则每次需要编译，比较耗时，（因为增加项会不断变化）。所以想通过在constant文件夹中以.H形式存储标量表达式，从而在turbulenceProperties中读取文件，以系数的形式修改表达式，由于系数是实时读取的，因此只需要编译一次。.H文件中就是一个标量表达式，命名为zetai.H
****问题：****如果仅有zetai.H，内容如下（I1 I2是已经在湍流模型中计算好了的体标量场）：

并在湍流模型中读取：

如何能将zetai中的内容转换为体标量场进行计算？

@李东岳 谢谢李老师的回复~ 教师节快乐~ 是一个标准的channel。但是雷诺数非常高，Re=80 million, 没有8000万网格，是NASA官网上的一个标准算例，独特的地方就是雷诺数非常大

设置了一个高雷诺数（Re=80million based on channel's height）的管道流动算例，在加密网格过程中，发现残差曲线出现了一些小幅度振荡，用的求解器是simpleFoam，
设置的simple算法收敛要求针对各个量是1e-4，所以1600迭代步就收敛了，但是相比于粗网格或者中等网格的结果，发现残差振荡更厉害，粗网格，中等网格与细网格残差对比如下：

而当使用可压缩求解器时，残差平缓很多，只是e的残差下降很慢。残差如下：

所以想请教大家，稳态不可压算例残差达到要求的标准是什么？除了各个量达到残差要求，有没有对曲线是振荡还是平缓有限制？

@李东岳 我重新定义了k_, 将其强制变成const，但是还是报错，错误如下：

error: no matching function for call to ‘ddt(const alphaField&, const rhoField&, const volScalarField&)’
                + fvc::ddt(alpha, rho, k_)
                  ~~~~~~~~^~~~~~~~~~~~~~~~

可见k_前面已经有const，但是还是不行，可见应该还是要修改alpha与rho的类型？但是我重新赋值定义alpha， rho为volScalarField, 会报初始化错误，不能直接这样改，修改方式如下：
const volScalarField alphat = this->alpha_; 再在fvc中调用fvc::ddt(alphat, rhot, k_)

@李东岳 谢谢李老师的回复~ 但是学生不太理解，volScalarField是继承于GeometricField的类型，k_是volScalarField类型，应该不冲突把？而且为什么fvm中也是这样定义的，也是传入这三个变量就可以用呢？

后续的故事一直在这里自问自答 把里面的量重新定义成带边界的量之后就可以用fvc了（除kSource()与fvOpetions(alpha, rho, k_)以外），现在还有一项没有搞定，请教下各位CFDer, 就是fvc::ddt(alpha, rho, k_)这一项，报错是

error: no matching function for call to ‘ddt(const alphaField&, const rhoField&, Foam::volScalarField&)’
                + fvc::ddt(alpha, rho, k_)
                  ~~~~~~~~^~~~~~~~~~~~~~~~

问题是fvc与fvm中此项定义是一样的，都是：

ddt
(
    const volScalarField& alpha,
    const volScalarField& rho,
    const GeometricField<Type, fvPatchField, volMesh>& vf
)
{
    return fv::ddtScheme<Type>::New
    (
        vf.mesh(),
        vf.mesh().ddtScheme
        (
            "ddt("
          + alpha.name() + ','
          + rho.name() + ','
          + vf.name() + ')'
        )
    ).ref().fvmDdt(alpha, rho, vf);
}

虽然说alphg的类型是const alphaField&, rho的类型是rhoField&, 它们俩是geometricOneField 不是 volScalarField，但是为什么fvm可以进行运算但是fvc不可以呢？

@fangyuanaza 我知道什么原因了，以tke方程第一项为例，在OpenFOAM V240版本中，fvc::SuSp(2.0/3.0*rho_*divU, k_) 这里面的量都是volScalarField, 定义rho_与divU的时候也是这样，但是在新版本V7中fvc::SuSp(2.0/3.0*alpha()*rho()*divU, k_) 中,这些量都是内部场的量，没有包含边界数据，改成fvc::SuSp(2.0/3.0*alpha*rho*divU, k_) ，并且divU的定义也从volScalarField::Internel 改成volScalarField, 就可以编译成功。但是这样太麻烦了吧，比如epsilonByk的传入变量也都是volScalarField::Internel形式的，要全部重新定义一边，这样感觉改起来太笨了，请教大家有什么好的改进方式么？从volScalarField::Internal的量，除了重新定义的方法，怎么能改成volScalarField，让显式离散fvc重新能用上，求tke方程残差呢？

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

@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 理解了，因此我应该取uTau的边界网格进行遍历

    forAll(uTau.boundaryFieldRef(),patchi)

或者

    forAll(uTau.boundaryField(),patchi)

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

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