OpenFOAM 中的壁面函数(二)

Giskard 60 0

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

1. 湍流动能 k 的壁面函数

OpenFOAM 中提供了两种k的壁面函数, kqRWallFunctionkLowReWallFunction

  • 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_wypl 的计算是一个迭代过程
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^2f

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

  • $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^2f 应该跟 \varepsilon\omega 那样(见后文),计算第一层网格内的值,并且考虑一个网格有多个边界面的情形。OpenFOAM 目前计算的是每一个边界面元上的值,不知道这两种方式对结果有多大影响。

发表评论 取消回复
表情 图片 链接 代码

分享