原有的湍流模型加上非线性项雷诺应力的问题

SHUKK

各位有经验的老师们，我想问一下加入修正雷诺应力项的修正湍流模型的问题。目前，我已经尝试了模仿ShihQuadraticKE湍流模型写法，重新写了一个新的湍流模型，具体修改的地方如下：

// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
void kEpsilonTBNN::correctNut()
{
    correctNonlinearStress(fvc::grad(U_));
}
void kEpsilonTBNN::correctNonlinearStress(const volTensorField& gradU)
{
    timeScale_=k_/epsilon_;

    // Linear (nut)
    nut_ = g1_*k_*timeScale_;
    nut_.correctBoundaryConditions();

    // nolinear(tau_NL)
    volSymmTensorField S(timeScale_*symm(gradU));
    volTensorField W(timeScale_*skew(gradU));

    nonlinearStress_ =
           2*k_
           *(
               g2_ * twoSymm(S&W)
             + g3_ * dev(innerSqr(S))
             + g4_ * dev(symm(W&W))
             + g5_ * twoSymm(W&innerSqr(S))
             + g6_ * dev(twoSymm(W&W&S))
             + g7_ * twoSymm(W&S&W&W)
             + g8_ * twoSymm(S&W&innerSqr(S))
             + g9_ * dev(twoSymm(innerSqr(S)&W&W))
             + g10_ * twoSymm(W&innerSqr(S)&W&W)
           );
}

这是我主要的修正项，即Aij

其中，g1到g10是我通过机器学习出来的数据，

g（i）跟T（i）都是有关Sij跟Ri的公式。目前，按照这种方法计算是发散的，其主要的原因是g（i）跟T（i）是随每次迭代求解出的U，k，epsilon发生变化的，但是只通过读取0文件下的g（i）场只是一个不变的标量，从而导致求解发散，因此需要将学习出来的神经网络结构耦合到OpenFOAM的环境中。因此我打算修改simpleFoam，linearViscousStrees和新加的reynoldsNet神经网络结构。

SHUKK

目前，是重新编写了simpleFoam，可以是在计算迭代中同时生成Sij跟Rij数据，为后面的linearViscousStrees提供相应数据，同时修改了动量方程的有关雷诺应力项。

tmp<fvVectorMatrix> tUEqn
    (
        fvm::div(phi, U)
      + MRF.DDt(U)
      + turbulence->divDevRhoReff(U, S, R)
     ==
      fvOptions(U)

我参照了之前的有关文献，同时修改了incompressible文件下的，incompressibleTurbulenceModel.H。

        //- Return the source term for the momentum equation using NN surrogate
        virtual tmp<fvVectorMatrix> divDevReff
        (
            volVectorField& U,
            volTensorField& S,
            volTensorField& R
        ) { }

incompressible/IncompressibleTurbulenceModel文件夹下的，IncompressibleTurbulenceModel.C和.H文件，添加了以下内容

template<class TransportModel>
Foam::tmp<Foam::fvVectorMatrix>
Foam::IncompressibleTurbulenceModel<TransportModel>::divDevReff
(
    volVectorField& U,
    volTensorField& S,
    volTensorField& R
) const
{
    return divDevRhoReff(U, S, R);
}

        //- Return the source term for the momentum equation
        virtual tmp<fvVectorMatrix> divDevReff
        (
            volVectorField& U,
            volTensorField& S,
            volTensorField& R
        ) const;

并修改了linearViscousStrees.C和.H文件，添加了以下内容

template<class BasicTurbulenceModel> 
Foam::tmp<Foam::fvVectorMatrix>
Foam::linearViscousStress<BasicTurbulenceModel>::divDevRhoReff
(
    volVectorField& U,
    volTensorField& S,
    volTensorField& R
) const
{
NN结构：a_dd
}
    return 
    (
        - fvm::laplacian(this->alpha_*this->rho_*this->nu(), U)
        - fvc::div((this->alpha_*this->rho_*this->nu())*dev2(T(fvc::grad(U))) - this->a_dd)
    );

a_dd就是我要加入的非线性项，里面的内容主要是链接了reynoldsNet神经网络。

SHUKK

目前，通过这些步骤我有几个问题。
1、修改的linearViscousStrees和IncompressibleTurbulenceModel.C和.H文件，我也需要像simpleFoam求解器一样修改名称和USER位置来编译，避免与原有的OF内容冲突吗？或是有什么特殊的编译方法？
2、我看之前有关的帖子中，说只要修改linearViscousStrees.C和.H文件来达到修正效果，那这修正IncompressibleTurbulenceModel.C和.H文件的意义是什么呢？
有做过相关内容的老师们，解答一下吗？

李东岳

g（i）跟T（i）都是有关Sij跟Ri的公式。目前，按照这种方法计算是发散的，其主要的原因是g（i）跟T（i）是随每次迭代求解出的U，k，epsilon发生变化的，但是只通过读取0文件下的g（i）场只是一个不变的标量，从而导致求解发散

这看起来是代码写的不对，原理是不应该是这样。所以我觉得你努力的方向没走对，最后一个帖子去debug意义也不大，因为你根本不需要那样做。你这个东西我比较感兴趣，你如果觉得你们老师对“联合培养研究生”这种合作感兴趣的话，我可以直接协助处理你OpenFOAM问题的代码。

SHUKK

@李东岳多谢李老师对我目前所做的项目有个极高的评价，由于我已经是研三了，这是我课题的最后一部分，“联合培养”项目开展时间有点晚了，多谢李老师的好意。对于我目前所做的东西，感觉李老师理解稍有偏差。g（i），也就是g1到g10是一个有关gardU，k，epsilon的函数，目前无法写出其显式表达，因此打算用神经网络结构来表示，所以按之前步骤只放到0文件夹下面就是一个不变量，迭代一次就应该更新了，每一次OF求解迭代出来的g（i）应该是不一样的。之前我是想跟python耦合的，即求出一步就用学习完成的神经网络更新一次g（i），但是对于高网络数量的模型计算就十分缓慢，因此我才打算将在python的神经网络结构转换为c++能识别的神经网络结构，并重新编写linearViscousStrees和IncompressibleTurbulenceModel来进行OF耦合。目前，修改的代码写的差不多了，就是不知linearViscousStrees和IncompressibleTurbulenceModel的有效编译步骤该如何进行，还需请教一下李老师。

李东岳

也就是g1到g10是一个有关gardU，k，epsilon的函数，目前无法写出其显式表达，因此打算用神经网络结构来表示

你可以让求解器每次都去constant文件夹里面去读取。不需要动linearViscousStrees和IncompressibleTurbulenceModel。代码动的越多越麻烦。

SHUKK

@李东岳李老师，是说可以模仿在constant/turbulenceProperties中修改湍流模型的Cmu这类常数数值的步骤一样吗？把g（i）的c++神经网络写到turbulenceProperties进去，让求解器读取，不用修改linearViscousStrees和IncompressibleTurbulenceModel。

李东岳

是的 :-)

SHUKK

@李东岳好的，李老师。我先按您的方法试试。

xpqiu

@SHUKK
你可以参考一下这个项目：https://github.com/furstj/myTurbulenceModels

里面有非线性湍流模型的例子，比如这个 EARSM 模型，https://github.com/furstj/myTurbulenceModels/tree/main/turbulenceModels/RAS/EARSM

你可以在这个 EARSM 模型的基础上来修改，基本上只需要其中的 correctNonlinearStress 函数就可以了。

SHUKK

@xpqiu 非常感谢，您的建议。我之前也是修改kepsilon湍流模型的，跟EARSM 模型一样，加了correctNonlinearStress函数。之前算发散的时候，我从g1到g10一个一个添加进去计算，发现只带入g1时候才能计算，并且得出来的效果与kOmegaSST的结果相似的，但是只要带入g2就开始发散了。后来，看之前做的过程时候发现，g（i），也就是g1到g10是一个有关gardU，k，epsilon的函数，目前无法写出其显式表达，应该用外接神经网络结构来计算更新。那为什么只带g1能算呢，因为g1*T1是等于OF的原有线性项假设，相当于没有改变其假设，只是修改了Cmu的值而已。
EARSM 模型的非线性项中，我看它的非线性项写法与我的类似，但是它的bate（i）是有明确的显式公式的，是能随着求解器迭代来计算更新的，而我的是隐式公式，目前没有明确公式。

    volScalarField beta1 = -N*(2.0*sqr(N) - 7.0*IIW) / Q;
    volScalarField beta3 = -12.0 * IV / (N * Q);
    volScalarField beta4 = -2.0 * (sqr(N)  - 2.0*IIW) / Q;
    volScalarField beta6 = -6.0 * N / Q;
    volScalarField beta9 =  6.0 / Q;
    volScalarField Cmu = - 0.5 * (beta1 + IIW * beta6);

    this->nut_ = Cmu * this->k_ * tau;
    this->nut_.correctBoundaryConditions();

    this->nonlinearStress_ = this->k_ * symm(
        beta3 * ( (W & W) - (1.0/3.0) * IIW * I )
        + beta4 * ( (S & W) - (W & S) )
        + beta6 * ( (S & W & W) + (W & W & S) - IIW * S - (2.0/3.0) * IV * I)
        + beta9 * ( (W & S & W & W) - (W & W & S & W) )

现在我要做的内容就让我的g（i）跟它bate（i）一样跟随迭代，为了达到我这个效果，不能单修改湍流模型，还要外加上我的c++神经网络来更新g（i）。目前就是编译的一些内容和步骤有点疑惑，需要请教。

xpqiu

@SHUKK
你的 g(i) 是用训练好的模型来计算对吧，那么其实你只需要用C++把你计算g(i)相关的代码写出来，写成湍流模型的成员函数，然后再在correctNonLinearStress里面调用这些函数就可以了。虽然g(i)不能写出显式的公式，但是也是代入一些流场变量，然后输出g(i)的数值吧，这样肯定也是可以写出可以调用的函数的。

SHUKK

@xpqiu 是的，应该是这个样的。但是有个问题，就是C++神经网络需要使用类似pytorch这样的环境，单加一个网络结构没办法识别计算，有那些类似pytorch这种的c++的神经网络编程语言吗？或者有什么神经网络语言是跟OF耦合的较好的呢？

xpqiu

@SHUKK
你需要 libtorch，这个是C++库，提供了跟 pytorch 一样的功能。也就是说，用pytorch训练的模型，可以做到利用 libtorch 提供的 API 读入到任何 C++ 程序中来用。
针对Open FOAM，这个项目 https://gitlab.com/tmaric/openfoam-ml 可以作为入门参考。

李东岳

是的，应该是这个样的。但是有个问题，就是C++神经网络需要使用类似pytorch这样的环境，单加一个网络结构没办法识别计算，有那些类似pytorch这种的c++的神经网络编程语言吗？或者有什么神经网络语言是跟OF耦合的较好的呢？

一楼不是都已经开始计算，并且发散了，最后怎么回归到如何进行计算的问题了？看起来好像还没跑起来。

你现在不是已经可以在每个时间步获取到g(i)的值了？目前仅仅需要把每个时间步的g(i)灌输到openfoam湍流模型里面就可以了？

SHUKK

@李东岳李老师，现在是可以在pytorch中每一步更新得到g（i），目前的问题是如何将pytorch中的神经网络能被OF读取识别，实现与OF耦合的过程。发散的是之前只在0文件下提供一个不变的标量场的方法。目前打算是将pytorch转变成c++的能识别神经网络库，g（i）的值可以通过这个库直接灌输到OF的每一个迭代步中。

SHUKK

@xpqiu 好的，多谢老师提供的资料。我之前打算是用caffe2来做为耦合的库，但是发现过于久远，可能有一些问题。我先了解一下libtorch的具体情况。

李东岳

我没深入到机器学习。从你这还学习了不少东西。按照我的理解，你通过机器学习+CFD的方法大体的路子就是：

每个步长首先通过pytorch获取g(i)，g(i)是一些列的标量
把g(i)的数值传输到OpenFOAM的非线性湍流模型里
每个时间步封闭，速度压力耦合求解
进入到下一个时间步，回归到第一步

是不是这个路子？如果是这样的话，为什么要将pytorch中的神经网络能被OF读取识别？

另外，机器学习是每个时间步都要更新g(i)么？？这是不是会特别慢。我之前浅显的看过一点，以为是一套计算域一套g(i)

SHUKK

@李东岳李老师理解的没问题。为什么要将pytorch中的神经网络能被OF读取识别？主要是为了解决计算速度过慢的问题。可以理解为通过机器学习，获得训练完成的神经网络，也就是一个有关g（i）隐式的公式，若能被OF识别并有效输出，那我就不用每次迭代计算调用pytorch来进行计算，在两个程序中来回计算，可以直接在OF求解中同步计算，相当于直接提供显式公式的效果。对于速度慢的问题，其实不用担心，我只需调用训练完成的神经网络结构即可，相当于用一个已知的公式获得我所需要的g（i）数据。

李东岳

@李东岳在原有的湍流模型加上非线性项雷诺应力的问题中说：

每个步长首先通过pytorch获取g(i)，g(i)是一些列的标量
把g(i)的数值传输到OpenFOAM的非线性湍流模型里
每个时间步封闭，速度压力耦合求解
进入到下一个时间步，回归到第一步

所以这是机器学习+CFD的一种方法。但是因为每一步都要训练出来一个g(i)，所以速度慢。

一个有关g（i）隐式的公式，若能被OF识别并有效输出，那我就不用每次迭代计算调用pytorch来进行计算

你们打算找出来g(i)的一个公式，写成$g(i)=f(\bfU)$类似的函数，然后流程就变成了：

最开始训练出$g(i)=f(\bfU)$函数表达式
进行CFD求解（在这一步里面，不需要花时间训练找参数）

我这样理解对么？

CFD中文网

原有的湍流模型加上非线性项雷诺应力的问题