CFD中文网

用你1楼提供的算例运行，发现在 subsetMesh这一步之后，得到的 frame 边界的 type 是 empty，并且 0 文件里面 U，p 场的 frame 边界的边界条件也变成了 empty。所以，之后的计算你这个 frame 边界都是以 empty 边界条件在运行，结果肯定就完全错了。
解决办法：在 subsetMesh 这一步之后把0文件里面那些场的边界条件修改过来，改成正确的设置再重新计算。
跑了几步，结果如下：

@旺财R @HIT_lly @youhaoyu @1064168551 已发

@bit_hypersonic
cellLimited 这个东西是一种 gradient limiter，你可以搜一下这个关键字了解一下其物理含义。也可以参考如下论文看看：Michalak, K., & Ollivier-Gooch, C. (2008). Limiters for unstructured higher-order accurate solutions of the euler equations. 46th AIAA Aerospace Sciences Meeting and Exhibit, January. https://doi.org/10.2514/6.2008-776

@Joker

1. 会使用设置的厚度，也就是说，二维算例 z 方向设置不同的厚度，虽然对流场没影响，但是会影响 force 算出来的力的值
2. 力的单位是 N，OpenFOAM 里面默认所有的单位都是 SI 量纲组成的单位。

@东方白杨 有一个 functionObject 叫 columnAverage，可以实现对一个方向进行平均。

@z597288 3楼用到的几个命令组合一下应该可以实现

可以用 LaTex 的语法来实现希腊字母输入。

k 和 epsilon 可以也设置为 slip，nut 设置为 calculated。

切完之后再用 clip ，或者 先用 Select Cells Through，选定你想显示的网格，再用 extract selection 把选定的网格提取出来单独显示即可。

1. 你的背景网格三个方向的尺寸差异太大了，z方向的尺寸只有其他两个方向的大概 1/6，这不利于 snap 。建议将背景网格尺寸调成一样，或者至少把差异缩小试试。
2. snap 阶段建议增加 featureSnap 的次数
3. meshQuality里面，maxNonOrtho 你设置的是45，建议调大一些，比如调到65。

试了一下，不完美，但是有改善。

这段代码跟使用SA湍流模型计算有一个差异是：这段代码是先求解 nuTilda 方程得到 nut，然后再用更新的nut来构建和求解 UEqn，而直接使用湍流模型的时候是先求解 UEqn，pEqn，然后再用修正的速度来构建和求解nuTildaEqn

@myheart OpenFOAM 里面求解理想气体（ideal gas，无粘性）流动的求解器是 rhoCentralFoam，你看这个求解器的算例壁面都用的是某种滑移边界。

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

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

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

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

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

@李东岳
只是知道用了这个函数就会引起计算错误，但是不确定是不是我其他代码有问题。
希望有大佬路过进来看看这个 bound 函数的实现到底为什么要这样。

@李东岳

不知道为什么这个 bound 函数里面要用到 average 这样的操作。
之所以关心这个问题，是因为这个 average 给我导致一个严重的bug

很多有界变量在计算过程中为了防止数值问题，需要限定下限值，在OpenFOAM中一般都是使用bound 函数，代码见：src/finiteVolume/cfdTools/general/bound/bound.C。然而这个 bound 函数的实现却有些让人摸不着头脑：

Foam::volScalarField&
Foam::bound(volScalarField& vsf, const dimensionedScalar& lowerBound)
{
    const scalar minVsf = min(vsf).value();

    if (minVsf < lowerBound.value())
    {
        Info<< "bounding " << vsf.name()
            << ", min: " << minVsf
            << " max: " << max(vsf).value()
            << " average: " << gAverage(vsf.primitiveField())
            << endl;

        vsf.primitiveFieldRef() = max
        (
            max
            (
                vsf.primitiveField(),
                fvc::average(max(vsf, lowerBound))().primitiveField()
              * pos0(-vsf.primitiveField())
            ),
            lowerBound.value()
        );

        vsf.boundaryFieldRef() = max(vsf.boundaryField(), lowerBound.value());
    }

    return vsf;
}

一般的想法是，直接把小于 lowerBound 的重置为 lowerBound 就行了。但是 OpenFOAM 的实现是，对于某个变量v，当 0 < v < lowerBound的，令 v=lowBound (与预期一致)；当 v < 0 时，pos0(-vsf.primitiveField()) 为 1，也就意味着用如下公式重新计算了下限值

fvc::average(max(vsf, lowerBound))().primitiveField()

这种做法有什么数学上的解释吗？

@李东岳
参考文献：
Knacke, T. (2013).
Potential effects of Rhie & Chow type interpolations in airframe
noise simulations. In: Schram, C., Dénos, R., Lecomte E. (ed):
Accurate and efficient aeroacoustic prediction approaches for
airframe noise, VKI LS 2013-03.

@Chen_hao
当这个模型用于多相流时，IOobject::groupName("omega", alphaRhoPhi.group()) 会变成类似 omega.air 这样的名字。