关于修改湍流模型中非线性雷诺应力项的问题

[Chen_hao]

@李东岳 谢谢老师对我问题的指导，理解了针对“动量方程”雷诺应力项的修改就是通过修改linearViscousStress.C。
我先介绍下我的研究方向：模拟三层液态金属组成的电池的不稳定性，随着电压增加时三层流体会发生不稳定性即湍流，目前组里采用DNS计算消耗计算资源太大，所以我需要采用湍流模型来尝试一下这个发生过程，我采用了基于磁流体修改后的K-omega模型，其需要加入磁动能源项与各项异性的雷诺应力，以下是我需要修改的非线性雷诺应力项：

我的问题是：雷诺应力也包含在kOmega.H中如下：

tmp<volTensorField> tgradU = fvc::grad(U);
    volScalarField G  
    (
        this->GName(),
        nut*(tgradU() && dev(twoSymm(tgradU())))
    );
    tgradU.clear()；

 tmp<fvScalarMatrix> omegaEqn	
    (
        fvm::ddt(alpha, rho, omega_) 
      + fvm::div(alphaRhoPhi, omega_)
      - fvm::laplacian(alpha*rho*DomegaEff(), omega_)
     ==
        gamma_*alpha*rho*G*omega_/k_
      - fvm::SuSp(((2.0/3.0)*gamma_)*alpha*rho*divU, omega_)
      - fvm::Sp(beta_*alpha*rho*omega_, omega_)
      + fvOptions(alpha, rho, omega_)  
      + fvm::SuSp(Comg_*S/k_, omega_) 
    );

其中G中包含线性雷诺应力，针对这个非线性雷诺应力是不是我需要同时：
1、修改动量方程中的，即通过修改linearViscousStress.C；
2、修改kOmega模型中的，即也要修改上面的G。

初学openfoam见识浅薄，希望各位老师和同学能不吝赐教。

[李东岳]

可以按照这个形式，把你需要植入的东西高亮一下么

[Chen_hao]

@李东岳

[李东岳]

你可以把公式翻译成我那种矢量的形式么 这是在openfoam里面做植入的第一步

[李东岳]

尤其是那个$S_k^{Diss}$，这个还看不出来能怎么写。不过要是这么植入的话，也可以，就是感觉有点硬植入的风格

[李东岳]

比如$\bfS=\nabla\bfU+\nabla\bfU^T$，$S_{ik}S_{ki}=\bfS\cdot\bfS$，其他的你得自己写一下，比如这个$\boldsymbol{\tau}$怎么用$\bfS, \bfU$等表示出来

[Chen_hao]

@李东岳
十分感谢老师的耐心回答！回复不及时请见谅，明白老师的意思了，方程转换好了（在编辑中稍等一会给您发过来），就您说的源项部分我单纯定义了一个标量场（语言基础差点），自己在尝试，请老师见谅。

[李东岳]

@chen_hao 嗯，你这种就是我说的类似硬植入，看起来有点那什么。如果Sdiff能转化就不需要这么植入了，如果不能转化那就这样了。就是看起来这种植入一般来说不多。

[Chen_hao]

@李东岳

关于雷诺应力张量形式，上面有错误，正确版如下：

不知道有没有错误，我也在写其代码，跟老师讨论。

[李东岳]

@chen_hao 你这个公式写的很好！我对比了一些kOmega。你应该把原本kOmega中的$G$变成你现在的这个$G=\tau:\bfU$就可以了。这个$\tau$通过张量公式里面相应的代码植入一下就好，也不需要硬植入，也就不需要修改linearViscousStress.C了。

[Chen_hao]

@李东岳
谢谢老师夸奖。您的意思是，只需要更换KOmega中的G，目的是求解更精准的$ν_t$，从而作用到动量方程的U求解，故不需要再修改动量方程中的雷诺应力项，就可以达到效果。不知道理解的对不对，我一直在思考既然修改的话，需要两个地方统一变换，才能达到效果，希望老师能解答这个迷惑。

[李东岳]

@chen_hao

只需要更换KOmega中的G，目的是求解更精准的ν, 从而作用到动量方程的U求解，故不需要再修改动量方程中的雷诺应力项，就可以达到效果

湍流模型的作用就是更改$\nu_t$这个值，所以你的理解是正确的，其他的不用改

[Chen_hao]

@李东岳
理解了，谢谢老师的耐心回答

[xpqiu]

@chen_hao
你可以看看 src/TurbulenceModels/incompressible/turbulentTransportModels/RAS/ShihQuadraticKE 这个湍流模型
这是一个非线性的湍流模型而且也只考虑了二阶非线性项，跟你要实现的那个类似。
这类模型除了修改 $\nu_t$ ，还需要把雷诺应力的非线性部分加上。

[李东岳]

对。

1. 你要是动kOmega模型，可以把$G$换成你的$\tau$，在linearStress.C那面还要处理一下devTau函数，要把非线性的贡献加上（$\tau$减去线性部分）。

2. 你要是动非线性类模型如ShihQuadraticKE，把$G$换成你的$\tau$。这样不需要特为处理devTau函数了。但要处理下omega方程（如果没有的话）。

[Chen_hao]

@xpqiu 谢谢老师的回答，ShihQuadraticKE模型是一个很好的参考！

这类模型除了修改 νt ，还需要把雷诺应力的非线性部分加上。

您所指的加上雷诺应力非线性部分就是UEqn.H中turbulence->divDevRhoReff(rho, U)项的修改对吗？

[xpqiu]

@chen_hao 在 关于修改湍流模型中非线性雷诺应力项的问题 中说：

您所指的加上雷诺应力非线性部分就是UEqn.H中turbulence->divDevRhoReff(rho, U)项的修改对吗？

不是，不需要修改这个函数，你仔细看看 ShihQuadraticKE 这个模型的继承关系，它继承的是 nonLinearEddyViscousity 类。这个类计算雷诺应力的方式不一样，除了线性项还有一个 nonLinearStress 项。你只要按照这个框架去写，把你的非线性项赋值给 nonLinearStress 就行了。

[Chen_hao]

@李东岳
我会采用老师提供的第一种思路进行修改，感激不尽~

[Chen_hao]

@xpqiu @李东岳

老师，改程序中有个关于语法的问题不太确定，ShiQuadraticKE.C 中程序没有采用模板化的，但是KOmaga.C是采用模板化的，在我更改KOmaga.H文件的继承类时，我上图这个写法对不对，直接在nonlinearEddyViscosityincompressible::RASModel 中类比eddyViscosity<RASModel<BasicTurbulenceModel>>。

[李东岳]

我觉得你不用这么写。你动继承肯定会有一些问题，目前也看不出来，不是很strightworrd。比较简单的是：

1. kOmega.H头文件不用动，你先去把tau植入到你的kOmega.C里面
2. 计算G
3. 去devTau函数里面把tau的非线性部分加上

你要是按照ShiQuadraticKE套路来，

1. 把k方程改成你的omega方程，不用动.H文件
2. 非线性项按照你的tau去写

第一种简单一些，但有一些硬植入。就是devTau函数这块
第二种更灵活。做完了就是非线性komega模型。

[Chen_hao]

@李东岳 @
谢谢老师指导，我最后按照方法二，完全按照他的格式来修改的，修改到最后我发现一个问题，如下

 public:

    //typedef typename BasicTurbulenceModel::alphaField alphaField;
    //typedef typename BasicTurbulenceModel::rhoField rhoField;
    //typedef typename BasicTurbulenceModel::transportModel transportModel;


    //- Runtime type information
    TypeName("MHDkOmega");


    // Constructors

        //- Construct from components
        MHDkOmega
        (
            //const alphaField& alpha,   ？？？？
           //const rhoField& rho,
            const geometricOneField& alpha,  
            const geometricOneField& rho,
            const volVectorField& U,
            const surfaceScalarField& alphaRhoPhi,
            const surfaceScalarField& phi,
            const transportModel& transport,
            const word& propertiesName = turbulenceModel::propertiesName,
            const word& type = typeName
        );

其中 alpha是定义为geometricOneField也就是定义为1，可我的模型是三层液体，采用了multiphaseInterFoam求解器又基于VOF，我才意识到是不是湍流模型一般不能适用于多相的情况，这个alpha能不能修改来适用于这个求解器，担心自己做了无用功，希望老师能帮忙解答

[李东岳]

从这个alpha来看，确实不能用于多相

[Chen_hao]

@李东岳

老师，像这种形式的alpha的定义，能用于多相的情况吗？

typedef typename BasicTurbulenceModel::alphaField alphaField
const alphaField& alpha

BasicTurbulenceModel这个通用类是代指的是KomegaSST？（就本例）

[Chen_hao]

@李东岳
请老师看看G这部分代码可不可以这样写，
矢量场Snl有点搞不清=和（）赋值，tgradU后面用不用加括号~

tmp<volTensorField> tgradU = fvc::grad(U);//tmp reload to Field
    volScalarField S2(2*magSqr(symm(tgradU())));
    
    tensor I(1,0,0,0,1,0,0,0,1);  
    volSymmTensorField S(twoSymm(tgradU));
    volTensorField W(twoSkew(tgradU));
    volTensorField Snl =
        1/omega_
       *(
           S&S-(1/3)*magSqr(W)*I
         + twoSymm(S&W)
         + (2.5)*(W&W-(1/3)*magSqr(W)*I)
        );
    
    volScalarField::Internal GbyNu(dev(twoSymm(tgradU()())+Snl) && tgradU()());
    volScalarField::Internal G(this->GName(), nut()*GbyNu);
    tgradU.clear();

[李东岳]

可以。我记得湍流方程右边都是只考虑网格点的作用，不需要考虑边界。因此你的G声明了内部场，没有问题。