CFD中文网

我也没看懂，楼主写得太复杂了，别人没功夫仔细看你的公式。我建议你自己用matlab算好每个mesh上的数值，然后写成csv格式用swak4Foam导入。

看上去像是ddt+div
查了一下代码:
https://github.com/OpenFOAM/OpenFOAM-3.0.x/blob/master/src/finiteVolume/finiteVolume/fvc/fvcDDt.C

template<class Type>
tmp<GeometricField<Type, fvPatchField, volMesh> >
DDt
(
    const surfaceScalarField& phi,
    const GeometricField<Type, fvPatchField, volMesh>& psi
)
{
    tmp<GeometricField<Type, fvPatchField, volMesh> > ddtDivPhiPsi
        = fvc::ddt(psi) + fvc::div(phi, psi);

    if (phi.mesh().moving())
    {
        return ddtDivPhiPsi - fvc::div(phi + phi.mesh().phi())*psi;
    }
    else
    {
        return ddtDivPhiPsi - fvc::div(phi)*psi;
    }
}

这是3.0的版本，不确定2.3是不是也是如此。从这个版本里看，就是ddt+div。

至于定义，应该是随体导数（Material derivative）。

我在知乎上写过一些gdb相关东西，你可以参考一下：
https://zhuanlan.zhihu.com/p/33117676

@vamoszhi

我一开始没看明白，是因为你这个公式里面的符号跟一般书上的不一样。

1. 看上去你这里的u'反而是average量，U_i才是湍动量。是么？
2. 公式7的雷诺应力tau是不是少了一项，2/3Kdelta_ij?
   
3. 你的意思是只在某个cell上套用公式9，其他cell继续用公式8么？
4. 我对能量方程不是很熟，公式8里面有P和p'，这就是两个变量了，你是不是还用了一些别的closure law来表达P变量了？

@wwzhao 谢谢。这个问题我现在已经搞清楚了，是为了增加矩阵的“对角占优”特性，提高收敛性。经过实验，假如修改成不加nu_t的公式，对于一些Co数较大的case就有可能计算发散。

LRR

tmp<fvVectorMatrix> LRR::divDevReff(volVectorField& U) const
{
    if (couplingFactor_.value() > 0.0)
    {
        return
        (
            fvc::div(R_ + couplingFactor_*nut_*fvc::grad(U), "div(R)")
          + fvc::laplacian
            (
                 (1.0 - couplingFactor_)*nut_,
                 U,
                 "laplacian(nuEff,U)"
            )
          - fvm::laplacian(nuEff(), U)
        );
    }
    else
    {
        return
        (
            fvc::div(R_)
          + fvc::laplacian(nut_, U, "laplacian(nuEff,U)")
          - fvm::laplacian(nuEff(), U)
        );
    }
}

tmp<fvVectorMatrix> UEqn
(
    fvm::ddt(U)
  + fvm::div(phi, U)
  + turbulence->divDevReff(U)
 ==
    fvOptions(U)
);


tmp<fvVectorMatrix> UEqn
(
    fvm::ddt(U)
  + fvm::div(phi, U)
  + turbulence->divDevReff(U)
 ==
    fvOptions(U)
);

且默认couplingFactor=0.

    //- Return the effective viscosity
        virtual tmp<volScalarField> nuEff() const
        {
            return tmp<volScalarField>
            (
                new volScalarField("nuEff", nut() + nu())
            );

@wwzhao LRR 求解器+pimpleFoam，我把两边的公式加起来得到的这个完整的公式。

其中R代表雷诺应力，S代表源项

我发现OpenFOAM2.4中把雷诺应力模型定义成了这样：

我不太理解为什么要这么写，以我的理解应该是这样：

请问有人能解开我的疑惑吗？

@李东岳 是的，我用python写了个回归函数，根据DNS给出的数据做回归，然后想往OpenFOAM里带。

@李东岳 多谢！我可以问问turbulence->correct()这句是干什么的吗？

@wwzhao 啊，sorry，没看清题目。

我想用python生成一个湍流模型的替代模型，放到OpenFOAM中，让OpenFOAM用我的python模型算动量方程中的雷诺应力项。请大家给我一些建议：
我打算修改pisoFoam.C

        fvVectorMatrix UEqn
        (
            fvm::ddt(U)
          + fvm::div(phi, U)
          + turbulence->divDevReff(U)
        );

改成

        fvVectorMatrix UEqn
        (
            fvm::ddt(U)
          + fvm::div(phi, U)
          - fvm::laplacian(nu+mynut, U)
        );

mynut是我算出来的一个tensorField，跟LRR里面一样处理。

但是我不确定是不是只改这里就行了，后面有句：

    turbulence->correct();

1. 我不知道这句有什么用，按照我的需求，是不是要把这句删掉？
2. 还有什么其他的地方需要改吗？

用pimpleFoam的话可以直接在fvoptions里添加momentum source项。

@xpqiu 哇，邱老师居然回答了我的问题。时间过去半年了，这个问题我已经搞懂了。我是群里的朱掌柜，前几周还问过你问题。

发现一个小问题：
http://dyfluid.com/icoFoam.html
对比公式1和公式5，似乎缺少了压力项。