OpenFOAM 中的壁面函数（二）

这是转载的文章，不是个人原创，,最初的作者我也不知道是谁了，就是觉得作者写的很好，我也想这么厉害[aru_9]
来看看计算湍动能 k 的壁面函数。

1. 湍流动能 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 目前计算的是每一个边界面元上的值，不知道这两种方式对结果有多大影响。