interFoam耦合相间的能量方程

[李东岳]

你先简单的弄一个纯对流方程试一下ddt(T)+div(phi,T)

[upc_ngh]

@李东岳 东岳老师好，我这个昨晚发现问题了，我原来的TEqn如下：

//计算错误
fvScalarMatrix TEqn
              (
                  fvm::ddt(rho*cp,Tcell)
                 +fvm::div(rhoPhiT,Tcell)
                 -fvm::laplacian(k,Tcell)
              );
             TEqn.solve(); 
//计算正确
fvScalarMatrix TEqn
              (
                  fvm::ddt(Tcell)
                 +fvm::div(phi,Tcell)
                 -fvm::laplacian(k/cp/rho,Tcell)
              );
             TEqn.solve(); 
//基金会版本
fvScalarMatrix TEqn
    (
        fvm::ddt(rhoCp, T)
      + fvm::div(rhoCpPhi, T, "div(phi,T)")
      - fvm::Sp(fvc::ddt(rhoCp) + fvc::div(rhoCpPhi), T)
      - fvm::laplacian(kappaEff, T,  "laplacian(kappa,T)")
      ==
        fluid.heatTransfer(T)
      + radiation->ST(T)
      + fvOptions(rhoCp, T)

上面求解是有问题的，后面我把各项的rho和cp同时约去，变成上面那样,这样就可以计算了，结果也是正常的，如下图。此外，我还有一个问题，如果想要考虑气液界面处的传热，是不是要在方程后面加上一项，就像基金会版本求解器icoReactingMultiphaseInterFoam中的代码中的 fluid.heatTransfer(T)，但具体这部分怎么做的，目前感觉没有思路也没找到相关的文献资料，还请东岳老师和各位大佬多多指点，谢谢！

[李东岳]

你为什么不直接用带传热的VOF，而是要自己植入呢

[upc_ngh]

@李东岳 非常感谢东岳老师的回复，接下来是（SOS，紧急求救！！！）首先说下模拟背景（OF9）：这个算例模拟辐射和对流加热下，熔渣受热熔化向下滴落这个过程如下图Alpha 0时刻 到 Alpha 1时刻 这个过程，觉得传热方程不复杂自己就植入进去了，以源项的形式考虑了熔渣和空气接触面处的传热，目前计算结果也没问题如图T 0时刻 到 T 1时刻，以上是考虑的第一部分。

第二部分考虑随着温度升高，熔渣的运动粘度不断变小，流动性变好，所以想植入nu与T之间的关系，但interFoam中nu是一个常数，根据代码发现，UEqn.H（第一段代码）中turbulence->divDevTau(rho, U)，使用到了nu，因此想在src库中找到关于divDevTau(rho, U)，把场变量Nu放进来，以为很简单，但接下来就开始了无穷尽的折磨。
（1）根据divDevTau()定义找到了位于src/MomentumTransportModels/incompressible/IncompressibleMomentumTransportModel/的函数定义（第二段代码），发现函数定义在反复横跳，就是不给表达式；
（2）根据网上资料发现在linearViscousStress.C文件下有关于****divDevTau(rho, U)****的定义，非常激动的准备把this->nuEff()，这个关于nu的函数直接改成Nu这个随温度变化的标量场（这里认为this->nuEff（）就是nu的场），结果毫无作用 ；
（3）既然上面没作用就想着直接把这段代码放在UEqn.H内（）如下UEqn.H的代码，发现量纲是对的，又很激动，结果能算，但是完全错误，可见上图的UY的对比，熔渣向下流动，结果我改完代码，速度是向上的。
（4）后面我把我采用笨办法把MomentumTransportModels中的所有代码都看了一遍，很不幸还是没找到关于这个歌函数的定义。。。。。。
（5）被逼疯的状态，国庆加班两天了，又到24点了，求求东岳老师和各位大佬拯救我吧，头发快被我抓掉光了，目前就是想找到 divDevTau(rho, U) ，这个SB函数的具体定义式，然后如何把nu这个场变量放进去，或者说this->nuEff()具体是啥，这我也找不到具体定义，都是虚函数，没有说明具体定义，唉要被逼疯了

    fvVectorMatrix UEqn
    (
        fvm::ddt(rho, U) + fvm::div(rhoPhi, U)
      + MRF.DDt(rho, U)
      + turbulence->divDevTau(rho, U)
      //- fvc::div((rho*NU)*dev2(T(fvc::grad(U))))%%（3）暴力放进来
    //  - fvm::laplacian(rho*NU, U)
     ==
       phaseChange.SU(rho, rhoPhi, U)
     + fvModels.source(rho, U)
    );

template<class TransportModel>
Foam::tmp<Foam::fvVectorMatrix>
Foam::IncompressibleMomentumTransportModel<TransportModel>::
divDevTau
(
    const volScalarField& rho,
    volVectorField& U
) const
{
    NotImplemented;
    
    return divDevSigma(U);
}

template<class TransportModel>
Foam::tmp<Foam::fvVectorMatrix>
Foam::IncompressibleMomentumTransportModel<TransportModel>::divDevSigma
(
    volVectorField& U
) const
{
    return divDevTau(U);
}

template<class TransportModel>
Foam::tmp<Foam::fvVectorMatrix>
Foam::IncompressibleMomentumTransportModel<TransportModel>::
divDevTau
(
    volVectorField& U
) const
{
    NotImplemented;
    return divDevSigma(U);
}

template<class BasicMomentumTransportModel>
Foam::tmp<Foam::fvVectorMatrix>
Foam::linearViscousStress<BasicMomentumTransportModel>::divDevTau
(
    const volScalarField& rho,
    volVectorField& U
) const
{

  volScalarField Nu = rho.db().objectRegistry::lookupObject<volScalarField>("NU");
    return
    (
      - fvc::div((this->alpha_*rho*this->nuEff()（改成Nu）)*dev2(T(fvc::grad(U))))
      - fvm::laplacian(this->alpha_*rho*this->nuEff()（改成Nu）, U)
    );
}

[李东岳]

@upc_ngh 在 interFoam耦合相间的能量方程 中说：

所以想植入nu与T之间的关系

你去viscosity模型那面，参考crossPowerLaw模型，植入一个自己的模型即可。不要管湍流模型那面的nueff()

[upc_ngh]

@李东岳 非常感谢李老师回复，现在最关键的问题是把nu与T之间的函数建立之后，无法将nu传递进divDevTau(rho, U)这个函数内部，因为根本找不到任何关于它的实质定义。这个问题有个前辈之前也发帖询问过：https://www.cfd-china.com/topic/3588/devreff-通过devrhoreff-实现其功能-rho-1-of4-x @Yongbo
（1）interFoam求解器内divDevTau(rho, U)指向的是 /IncompressibleMomentumTransportModel.C下的divDevTau(rho, U)这个函数，但是根本没找到任何相关的divDevTau(rho, U)的具体定义；
（2）反而在linearViscousStress.C下存在divDevTau(rho, U)，如下，但发现改变这个没任何反应，说明调用的也不是这个，后面我把IncompressibleMomentumTransportModel这个类下的所有有关divDevTau这个函数都改了一遍，发现都没任何反应，人完全懵了，就想知道interFoam中divDevTau(rho, U)具体是怎么实现的，虽然在IncompressibleMomentumTransportModel下，但是根本没有任何定义，人都被代码搞麻了

template<class BasicMomentumTransportModel>
Foam::tmp<Foam::fvVectorMatrix>
Foam::linearViscousStress<BasicMomentumTransportModel>::divDevTauBaby
(
    const volScalarField& rho,
    volVectorField& U
) 
{

  volScalarField Nu = rho.db().objectRegistry::lookupObject<volScalarField>("NU");%%这是我写的想要将变化的Nu传递进来
    return
    (
      - fvc::div((this->alpha_*rho*Nu)*dev2(T(fvc::grad(U))))
      - fvm::laplacian(this->alpha_*rho*Nu, U)
    );
}

[李东岳]

@upc_ngh 在 interFoam耦合相间的能量方程 中说：

现在最关键的问题是把nu与T之间的函数建立之后，无法将nu传递进divDevTau(rho, U)这个函数内部

为什么，你的意思所有的粘度模型visosity再湍流下都是失效的么

[upc_ngh]

@李东岳 目前应该是还没找对divDevTau(rho, U)具体位置，不知道把nu传给谁，网上说传给 linearViscousStress.C，但没任何反应，我昨天有个笨办法，在OF9下的MomentumTransportModels这个文件夹下关于divDevTau()所有定义内都加了Info输出某个量，用来监测interFoam中的divDevTau()到底调用了哪个函数，如下，但是没发现任何调用的迹象，我感觉很懵，明明定义在这个文件夹下，咋会找不到它的具体调用 ，东岳老师Help !!!现在完全没思路了

%%这是非线性的divDevTau函数定义在  nonlinearEddyViscosity.C内
template<class BasicMomentumTransportModel>
Foam::tmp<Foam::fvVectorMatrix>
Foam::nonlinearEddyViscosity<BasicMomentumTransportModel>::divDevTau
(
    const volScalarField& rho,
    volVectorField& U
) const
{

Info<< " 13 "<<endl;
    return
    (
        fvc::div(rho*nonlinearStress_)
      + eddyViscosity<BasicMomentumTransportModel>::divDevTau(rho, U)
    );
}

%%这是非线性的divDevTau函数定义在  linearViscousStress.C内
template<class BasicMomentumTransportModel>
Foam::tmp<Foam::fvVectorMatrix>
Foam::linearViscousStress<BasicMomentumTransportModel>::divDevTau
(
    const volScalarField& rho,
    volVectorField& U
) const
{

  //volScalarField Nu = rho.db().objectRegistry::lookupObject<volScalarField>("NU");
  Info<<"8"<<endl;
    return
    (
      - fvc::div((this->alpha_*rho*this->nuEff())*dev2(T(fvc::grad(U))))
      - fvm::laplacian(this->alpha_*rho*this->nuEff(), U)
    );
}

[李东岳]

1. 编译完了lib之后再编译一下求解器确保正确

2. 你不需要管湍流模型的调用，你要改的是viscosity模型而不是湍流模型

[upc_ngh]

@李东岳 非常感谢李老师不辞辛苦的点拨，现在有了一点小进展，为避免之后的同学犯错总结一下：
（1）之前错误点在于：通过lookup寻找随温度变化的nu场，再将 linearViscousStress.C中的divDevTau()下的this->nuEff()转换为nu场如下代码，这时仅仅编译 linearViscousStress.C上层的momentumTransportModels是不够的，还需要编译incompressible中的divDevTau()函数（尽管并没有修改这个函数），再编译求解器才可以。证明了interFoam中的divDevTau()中的函数在层流下执行的就是 linearViscousStress.C中的divDevTau()。
（2）但现在又存在一个问题，按照东岳老师的interFoam的解析并对照下面的代码可以看出， this->nuEff()是等价与nu即图中的v(运动粘度)

但是如果简单的将this->nuEff()代替成nu是不对的，计算出来的流场是错误的（熔渣向下滴落结果速度是向上的）如下图，其中nu的值都是根据气液相分数乘以分别的nu值得到，不知道哪里错了，或许这里的this->nuEff并不是简单的nu的数值，而是对其进行了处理，但也没找到相关定义，真是举步维艰，心累唉

%  linearViscousStress.C中的divDevTau()函数返回值
volScalarField Nu = rho.db().objectRegistry::lookupObject<volScalarField>("nu");
    return
    (
      - fvc::div((this->alpha_*rho*this->nuEff())*dev2(T(fvc::grad(U))))
      - fvm::laplacian(this->alpha_*rho*this->nuEff(), U)
    );

%  修改后linearViscousStress.C中的divDevTau()
volScalarField Nu = rho.db().objectRegistry::lookupObject<volScalarField>("nu");
    return
    (
      - fvc::div((this->alpha_*rho*Nu)*dev2(T(fvc::grad(U))))
      - fvm::laplacian(this->alpha_*rho*Nu, U)
    );

[李东岳]

@upc_ngh 在 interFoam耦合相间的能量方程 中说：

this->nuEff并不是简单的nu的数值

这个是两个年度场的混合，类似密度的数值

[upc_ngh]

@李东岳 感谢李老师的回复，根据查阅的资料：this->nuEff() = nu() + nut() interFoam求解器使用的this->nuEff()应该是laminarModel下的Stokes.C这个模型，因为在求解器计算时输出器选用的stress model 是Stokes model，那divDevTau()函数下的this->nuEff()就如Stokes.C内所描述如下，就是指向nu()并不是为this->nut() + this->nu()，这里就很疑惑，这里将nu()换成自己设定的Nu(T)，就得不到合理的计算结果（流场分布）。此外，我选用的是层流在laminarModel.C内也将nut初始化为0,如下代码，所以不知道interFoam这里调用了哪里的的nuEff()，可能不是目前我认知的Stokes.C内部的

template<class BasicMomentumTransportModel>
tmp<volScalarField> Stokes<BasicMomentumTransportModel>::nuEff() const
{
    return volScalarField::New
    (
        IOobject::groupName("nuEff", this->alphaRhoPhi_.group()),
        this->nu()
    );
}

template<class BasicMomentumTransportModel>
Foam::tmp<Foam::volScalarField>
Foam::laminarModel<BasicMomentumTransportModel>::nut() const
{
    return volScalarField::New
    (
        IOobject::groupName("nut", this->alphaRhoPhi_.group()),
        this->mesh_,
        dimensionedScalar(dimViscosity, 0)
    );
}