2016-04-25

OpenFOAM 中的壁面函数（二）

这篇来看看计算湍动能 $k$ 的壁面函数。

2. 湍流动能 $k$ 的壁面函数

OpenFOAM 中提供了两种 $k$ 的壁面函数， kqRWallFunction 和 kLowReWallFunction 。

kqRWallFunction
其实就是 zeroGradient ，无需多言。除非使用 $v^2\text{-}f$ 模型，一般情况下 $k$ 应该使用这个边界条件。

kLowReWallFunction
这个壁面函数应该是可以用于低雷诺数模型的。该壁面函数继承自 fixedValue ：

class kLowReWallFunctionFvPatchScalarField
:
    public fixedValueFvPatchField<scalar>
{
protected:
       //- Cmu coefficient
        scalar Cmu_;

        //- Von Karman constant
        scalar kappa_;

        //- E coefficient
        scalar E_;

        //- Ceps2 coefficient
        scalar Ceps2_;

        //- Y+ at the edge of the laminar sublayer
        scalar yPlusLam_;
        ......
        
kLowReWallFunctionFvPatchScalarField::kLowReWallFunctionFvPatchScalarField
(
    const fvPatch& p,
    const DimensionedField<scalar, volMesh>& iF,
    const dictionary& dict
)
:
    fixedValueFvPatchField<scalar>(p, iF, dict),
    Cmu_(dict.lookupOrDefault<scalar>("Cmu", 0.09)),
    kappa_(dict.lookupOrDefault<scalar>("kappa", 0.41)),
    E_(dict.lookupOrDefault<scalar>("E", 9.8)),
    Ceps2_(dict.lookupOrDefault<scalar>("Ceps2", 1.9)),
    yPlusLam_(yPlusLam(kappa_, E_))
    {
        checkType();
    }

}

核心的函数是以下两个：

scalar kLowReWallFunctionFvPatchScalarField::yPlusLam
(
    const scalar kappa,
    const scalar E
)
{
    scalar ypl = 11.0;

    for (int i=0; i<10; i++)
    {
        ypl = log(max(E*ypl, 1))/kappa;
    }

    return ypl;
}
void kLowReWallFunctionFvPatchScalarField::updateCoeffs()
{
    if (updated())
    {
        return;
    }

    const label patchI = patch().index();

    const turbulenceModel& turbulence =
        db().lookupObject<turbulenceModel>("turbulenceModel");
    const scalarField& y = turbulence.y()[patchI];

    const tmp<volScalarField> tk = turbulence.k();
    const volScalarField& k = tk();

    const tmp<volScalarField> tnu = turbulence.nu();
    const scalarField& nuw = tnu().boundaryField()[patchI];

    const scalar Cmu25 = pow025(Cmu_);

    scalarField& kw = *this; // 更新 kw 相当于更新壁面上的 k 值。

    // Set k wall values
    forAll(kw, faceI)
    {
        label faceCellI = patch().faceCells()[faceI];

        scalar uTau = Cmu25*sqrt(k[faceCellI]);

        scalar yPlus = uTau*y[faceI]/nuw[faceI];

        if (yPlus > yPlusLam_)
        {
            scalar Ck = -0.416;
            scalar Bk = 8.366;
            kw[faceI] = Ck/kappa_*log(yPlus) + Bk;
        }
        else
        {
            scalar C = 11.0;
            scalar Cf = (1.0/sqr(yPlus + C) + 2.0*yPlus/pow3(C) - 1.0/sqr(C));
            kw[faceI] = 2400.0/sqr(Ceps2_)*Cf;
        }

        kw[faceI] *= sqr(uTau);
    }

    fixedValueFvPatchField<scalar>::updateCoeffs();

    // TODO: perform averaging for cells sharing more than one boundary face
}

先在函数里计算 ypl 的值， updateCoeffs 函数里根据 yPlus 与这个 ypl 的值来相对大小而采取不同的方法来计算壁面上的 $k_w$。 ypl 的计算是一个迭代过程
$$
ypl = \frac{\log(\max(E*ypl,1.0))}{\kappa}
$$
初始值为 ypl = 11.0，迭代10次，最终结果应该是 ypl = 11.5301073043272。
$y^+$ 定义为：
$$
u_\tau = C_\mu^{1/4 }\sqrt{k_c} \
y^+ = \frac{u_\tau \cdot y}{\nu_w}
$$
壁面上的k计算方法如下：如果 $y^+ > ypl$，则
$$
k^+ _w = \frac{C_k}{\kappa}\ln(y^+) + B_k
$$
否则
$$
k^+ _w = \frac{2400}{C_{eps2}^2}\cdot \left[ \frac{1}{(y^+ + C)^2} + \frac{2y^+}{C^3} - \frac{1}{C^2}\right ]
$$
最终，壁面上的值为 $k_w=k^+ _w u_\tau ^2 =k^+ _w C_\mu^{1/2}k_c$ 。
以上公式中，下标 $c$ 表示壁面单元所述网格的值，下标 $w$ 表示当前壁面上的值。
这个壁面函数参考文献 “Kalitzin, G., Medic, G., Iaccarino, G., Durbin, P., 2005. Near-wall behavior of RANS turbulence models and implications for wall functions. J. Comput. Phys. 204, 265–291. doi:10.1016/j.jcp.2004.10.018”，是为 $v^2\text{-}f$ 模型设计的。

$v^2$ 和 $f$ 的壁面函数

上面提到了 $v^2\text{-}f$ 模型，所以这里顺便来看看$v^2$ 和 $f$ 的壁面函数。这里参考的也是上面提到的那篇参考文献。

$v^2$ 的壁函数

forAll(v2, faceI)
    {
        label faceCellI = patch().faceCells()[faceI];

        scalar uTau = Cmu25*sqrt(k[faceCellI]);

        scalar yPlus = uTau*y[faceI]/nuw[faceI];

        if (yPlus > yPlusLam_)
        {
            scalar Cv2 = 0.193;
            scalar Bv2 = -0.94;
            v2[faceI] = Cv2/kappa_*log(yPlus) + Bv2;
        }
        else
        {
            scalar Cv2 = 0.193;
            v2[faceI] = Cv2*pow4(yPlus);
        }

        v2[faceI] *= sqr(uTau);
    }

    fixedValueFvPatchField<scalar>::updateCoeffs();

yPlus > yPlusLam_ 时，
$$
v^2 = u_\tau^2 \cdot \left[ \frac{C_{v2}}{\kappa}\ln(y^+) + B_{v2} \right]
$$
与文献中的无量纲形式 $(\overline{v^2})^{^+} = \frac{C_{v2}}{\kappa}\ln(y^+) + B_{v2} $ 一致。

yPlus < yPlusLam\_ 时，
$$
v^2 = u_\tau^2 \cdot C_{v2}(y^+)^2
$$
与无量纲形式 $(\overline{v^2})^{^+} = C_{v2}(y^+)^2$ 一致。

$f$ 的壁函数

forAll(f, faceI)
    {
        label faceCellI = patch().faceCells()[faceI];

        scalar uTau = Cmu25*sqrt(k[faceCellI]);

        scalar yPlus = uTau*y[faceI]/nuw[faceI];

        if (yPlus > yPlusLam_)
        {
            scalar N = 6.0;
            scalar v2c = v2[faceCellI];
            scalar epsc = epsilon[faceCellI];
            scalar kc = k[faceCellI];

            f[faceI] = N*v2c*epsc/(sqr(kc) + ROOTVSMALL);
            f[faceI] /= sqr(uTau) + ROOTVSMALL;
        }
        else
        {
            f[faceI] = 0.0;
        }
    }

yPlus > yPlusLam_ 时，
$$
f = \frac{N \cdot v^2\cdot \varepsilon}{k^2 u_\tau^2}
$$
这似乎与文献中的无量纲形式
$$
f^+ = N \frac{(\overline{v^2})^{^+}}{(k^+)^2}\varepsilon^+
$$
不一致！是 bug 还是我推导出错了？存疑…

yPlus < yPlusLam_ 时，文献给出的公式是
$$
f^+ = \frac{-4(6-N)(\overline{v^2})^{^+}}{\varepsilon^+ (y^+)^4}
$$
当 N=6 时，可以得到 $f^+ = 0$ 。

按理说，$v^2$ 和 $f$ 应该跟 $\varepsilon$ 和 $\omega$ 那样（见后文），计算第一层网格内的值，并且考虑一个网格有多个边界面的情形。OpenFOAM 目前计算的是每一个边界面元上的值，不知道这两种方式对结果有多大影响。

Giskard's CFD Learning Tricks

CFD and Scientific Computing

OpenFOAM 中的壁面函数（二）

2. 湍流动能 $k$ 的壁面函数

$v^2$ 和 $f$ 的壁面函数