bound 函数的实现方法

xpqiu

很多有界变量在计算过程中为了防止数值问题，需要限定下限值，在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()

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

李东岳

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

$v_{celli}=\frac{\sum (|\bfS_f| v_f)}{\sum |\bfS_f|}$

应该是基金会自己拍脑袋定的。

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

大佬这不就是average的实现么？

考虑这么几个网格|3|-1|2|

max(vsf, lowerBound)变成了|3|0|2|

fvc::average(max(vsf, lowerBound))().primitiveField()变成了|3|1.25|2|

xpqiu

@李东岳

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

李东岳

大佬那面玩的高端。你说是bug，那绝对是bug。跟openfoam.org提呗。

xpqiu

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

李东岳

bound 函数的实现到底为什么要这样。

这个看起来是个很简单的数值操作。不是从物流守恒定律搞出来的。这就是Henry自己搞的。他很多东西都自己搞。

average这个函数我最大的印象就是dynSmagorinsky模型，Henry搞这么个东西进去，整个场的动态Cs变成了一个标量，后来让alberto怼有bug，Henry直接给dynSmagorinsky模型给删了。到现在OpenFOAM也没有动态Smagorinsky模型。 https://bugs.openfoam.org/view.php?id=99

cyberk

@李东岳

帮李老师补充一下OpenFOAM里面average的实现，对于vol场要先插值成surface场

template<class Type>
tmp<GeometricField<Type, fvPatchField, volMesh>>
average
(
    const GeometricField<Type, fvPatchField, volMesh>& vtf
)
{
    return fvc::average(linearInterpolate(vtf));
}

所以|3|0|2|到average函数里面以后要先插值成

|  3  |  0  |  2  | 
     1.5    1

再基于

$v_{celli}=\frac{\sum (|\bfS_f| v_f)}{\sum |\bfS_f|}$
得到中间网格的值为1.25

CFD中文网

bound 函数的实现方法