2016-04-25

OpenFOAM 中的壁面函数（四）

这篇来看看可能是最关键的 $\nu_t$ 的壁面函数。

5. 湍流粘度 $\nu_t$ 的壁面函数

这个类型的壁面函数，结构比较简单，计算的是每一个壁面边界面上的湍流粘度 $\nu_t$。
nutWallFunction 是虚基类，其中定义了一个纯虚函数 calcNut

1	virtual tmp<scalarField> calcNut() const = 0;

并且在 updateCoeffs 函数中，将 calcNut 的返回值赋值给边界面

void nutWallFunctionFvPatchScalarField::updateCoeffs()
{
    if (updated())
    {
        return;
    }

    operator==(calcNut());

    fixedValueFvPatchScalarField::updateCoeffs();
}

这样，在具体的那些计算 $\nu_t$ 的壁面函数中，只需要看 calcNut 的返回值就可以了。

(1). nutkWallFunction

tmp<scalarField> nutkWallFunctionFvPatchScalarField::calcNut() const
{
    const label patchi = patch().index();
    const turbulenceModel& turbModel =
        db().lookupObject<turbulenceModel>("turbulenceModel");
    const scalarField& y = turbModel.y()[patchi];
    const tmp<volScalarField> tk = turbModel.k();
    const volScalarField& k = tk();
    const tmp<volScalarField> tnu = turbModel.nu();
    const volScalarField& nu = tnu();
    const scalarField& nuw = nu.boundaryField()[patchi];
    const scalar Cmu25 = pow025(Cmu_);
    tmp<scalarField> tnutw(new scalarField(patch().size(), 0.0));
    scalarField& nutw = tnutw();
    forAll(nutw, faceI)
    {
        label faceCellI = patch().faceCells()[faceI];

        scalar yPlus = Cmu25*y[faceI]*sqrt(k[faceCellI])/nuw[faceI];

        if (yPlus > yPlusLam_)
        {
            nutw[faceI] = nuw[faceI]*(yPlus*kappa_/log(E_*yPlus) - 1.0);
        }
    }
    return tnutw;
}

这里，仍然是分情况处理
yPlus < yPlusLam_ 时，壁面上的 nut 设为0；
yPlus > yPlusLam_ 时
$$
\nu_t = \nu \cdot \left( \frac{\kappa y^+}{\ln(Ey^+)}-1 \right)
$$
这里实现的其实就是标准壁面函数。理论上讲，这里的计算只在粘性底层和对数区是有效的，所以，使用这个壁面条件的时候，要尽量壁面网格落在过渡区，否则可能会引入较大误差。

顺带提一下，这里还定义了一个 yPlus 函数，用来计算 $y^+$，这个函数在这里没有调用，不过在其他代码中需要 $y^+$ 的时候会调用这个函数。比如，计算 $y^+$ 的应用 yPlusRAS 就是调用这里的 yPlus 函数来计算 $y^+$。

tmp<scalarField> nutkWallFunctionFvPatchScalarField::yPlus() const
{
    const label patchi = patch().index();
    const turbulenceModel& turbModel =
        db().lookupObject<turbulenceModel>("turbulenceModel");
    const scalarField& y = turbModel.y()[patchi];

    const tmp<volScalarField> tk = turbModel.k();
    const volScalarField& k = tk();
    tmp<scalarField> kwc = k.boundaryField()[patchi].patchInternalField();
    const tmp<volScalarField> tnu = turbModel.nu();
    const volScalarField& nu = tnu();
    const scalarField& nuw = nu.boundaryField()[patchi];

    return pow025(Cmu_)*y*sqrt(kwc)/nuw;
}

(2). nutUWallFunction

tmp<scalarField> nutUWallFunctionFvPatchScalarField::calcNut() const
{
    const label patchi = patch().index();
    const turbulenceModel& turbModel =
        db().lookupObject<turbulenceModel>("turbulenceModel");
    const fvPatchVectorField& Uw = turbModel.U().boundaryField()[patchi];
    const scalarField magUp(mag(Uw.patchInternalField() - Uw));
    const tmp<volScalarField> tnu = turbModel.nu();
    const volScalarField& nu = tnu();
    const scalarField& nuw = nu.boundaryField()[patchi];
    tmp<scalarField> tyPlus = calcYPlus(magUp);
    scalarField& yPlus = tyPlus();
    tmp<scalarField> tnutw(new scalarField(patch().size(), 0.0));
    scalarField& nutw = tnutw();

    forAll(yPlus, facei)
    {
        if (yPlus[facei] > yPlusLam_)
        {
            nutw[facei] =
                nuw[facei]*(yPlus[facei]*kappa_/log(E_*yPlus[facei]) - 1.0);
        }
    }
    return tnutw;
}

tmp<scalarField> nutUWallFunctionFvPatchScalarField::calcYPlus
(
    const scalarField& magUp
) const
{
    const label patchi = patch().index();
    const turbulenceModel& turbModel =
        db().lookupObject<turbulenceModel>("turbulenceModel");
    const scalarField& y = turbModel.y()[patchi];
    const tmp<volScalarField> tnu = turbModel.nu();
    const volScalarField& nu = tnu();
    const scalarField& nuw = nu.boundaryField()[patchi];
    tmp<scalarField> tyPlus(new scalarField(patch().size(), 0.0));
    scalarField& yPlus = tyPlus();
    forAll(yPlus, facei)
    {
        scalar kappaRe = kappa_*magUp[facei]*y[facei]/nuw[facei];
        scalar yp = yPlusLam_;
        scalar ryPlusLam = 1.0/yp;
        int iter = 0;
        scalar yPlusLast = 0.0;
        do
        {
            yPlusLast = yp;
            yp = (kappaRe + yp)/(1.0 + log(E_*yp));

        } while (mag(ryPlusLam*(yp - yPlusLast)) > 0.01 && ++iter < 10 );

        yPlus[facei] = max(0.0, yp);
    }
    return tyPlus;
}

这个壁函数的 $y^+$ 的计算方式跟 nutkWallFunction 有点区别。经过摸索，这里 calcYPlus 函数中的那段 do ... while 循环的原理如下：
对数律可以表达如下：
$$
U^+ = \frac{U_p}{u_\tau}=\frac{1}{\kappa}\ln(Ey^+)
$$
其中 $U_p$ 等于壁面上的速度减去壁面所属网格中心的速度。
经过简单变形
$$
\frac{U_p}{ y u_\tau/\nu }\cdot (y/\nu)=\frac{U_p}{ y^+}\cdot (y/\nu)=\frac{1}{\kappa}\ln(Ey^+)
$$
整理得
$$
y^+ \ln(Ey^+) - \frac{\kappa y U_p}{\nu}=0
$$
这是一个 $y^+$ 的一元方程，可以通过牛顿迭代来求解
$$
y^+_{n+1} = y^+_{n} - \frac{f(y^+)}{f^{\prime}(y+)} = y^+_{n}-\frac{y_n^+ \ln(Ey_n^+) - \frac{\kappa y U_p}{\nu}}{1+\ln(Ey_n^+)} = \frac{y_n^+ + \frac{\kappa y U_p}{\nu}}{1+\ln(Ey_n^+)}
$$
上面代码里的 do ... while 循环，正是在做这个迭代求解，初始值选择的是 yPlusLam，这个值在前面提过了。
求出 $y^+$ 以后，$\nu_t$ 计算如下
$$
\nu_t = \nu \cdot \left( \frac{\kappa y^+}{\ln(Ey^+)}-1 \right)
$$
与 nutkWallFunction 形式是一样的。

这个壁面函数，求壁面上的 $\nu_t$ 时使用的对数律方程，所以，理论上这个壁面函数应该只适用于第一层网格落在对数层的情形。

(3). nutLowReWallFunction
这个壁面函数直接将壁面上的 $\nu_t$ 的值设为0。

tmp<scalarField> nutLowReWallFunctionFvPatchScalarField::calcNut() const
{
    return tmp<scalarField>(new scalarField(patch().size(), 0.0));
}

$y^+$ 的计算也值得注意：

tmp<scalarField> nutLowReWallFunctionFvPatchScalarField::yPlus() const
{
    const label patchi = patch().index();
    const turbulenceModel& turbModel =
        db().lookupObject<turbulenceModel>("turbulenceModel");
    const scalarField& y = turbModel.y()[patchi];
    const tmp<volScalarField> tnu = turbModel.nu();
    const volScalarField& nu = tnu();
    const scalarField& nuw = nu.boundaryField()[patchi];
    const fvPatchVectorField& Uw = turbModel.U().boundaryField()[patchi];

    return y*sqrt(nuw*mag(Uw.snGrad()))/nuw;
}

$$
y^+ = \frac{y\sqrt{\nu \cdot |\frac{U_w-U_c}{d}|}}{\nu}
$$
注意由于 $\nu_t = 0$ ，所以 $\frac{\tau_w}{\rho} = \nu \cdot |\frac{U_w-U_c}{d}|$，所以，$\sqrt{\nu \cdot |\frac{U_w-U_c}{d}|}=\sqrt{\frac{\tau_w}{\rho}}=u_\tau$ 。

(4). nutUSpaldingWallFunction
这个壁函数基于 Spalding 提出的一个拟合的 $y^+$ 与 $u^+$ 的关系式，见文献 A Single Formula for the “Law of the Wall” 。
$$
y^+ = u^+ + \frac{1}{E}\left[ e^{\kappa u^+} -1-\kappa u^+ -\frac{1}{2}(\kappa u^+)^2 - \frac{1}{6}(\kappa u^+)^3 \right]
$$

tmp<scalarField> nutUSpaldingWallFunctionFvPatchScalarField::calcNut() const
{
    const label patchI = patch().index();

    const turbulenceModel& turbModel =
        db().lookupObject<turbulenceModel>("turbulenceModel");
    const fvPatchVectorField& Uw = turbModel.U().boundaryField()[patchI];
    const scalarField magGradU(mag(Uw.snGrad()));
    const tmp<volScalarField> tnu = turbModel.nu();
    const volScalarField& nu = tnu();
    const scalarField& nuw = nu.boundaryField()[patchI];

    return max
    (
        scalar(0),
        sqr(calcUTau(magGradU))/(magGradU + ROOTVSMALL) - nuw
    );
}

tmp<scalarField> nutUSpaldingWallFunctionFvPatchScalarField::calcUTau
(
    const scalarField& magGradU
) const
{
    const label patchI = patch().index();

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

    const fvPatchVectorField& Uw = turbModel.U().boundaryField()[patchI];
    const scalarField magUp(mag(Uw.patchInternalField() - Uw));

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

    const scalarField& nutw = *this;

    tmp<scalarField> tuTau(new scalarField(patch().size(), 0.0));
    scalarField& uTau = tuTau();

    forAll(uTau, faceI)
    {
        scalar ut = sqrt((nutw[faceI] + nuw[faceI])*magGradU[faceI]);

        if (ut > ROOTVSMALL)
        {
            int iter = 0;
            scalar err = GREAT;

            do
            {
                scalar kUu = min(kappa_*magUp[faceI]/ut, 50);
                scalar fkUu = exp(kUu) - 1 - kUu*(1 + 0.5*kUu);

                scalar f =
                    - ut*y[faceI]/nuw[faceI]
                    + magUp[faceI]/ut
                    + 1/E_*(fkUu - 1.0/6.0*kUu*sqr(kUu));

                scalar df =
                    y[faceI]/nuw[faceI]
                  + magUp[faceI]/sqr(ut)
                  + 1/E_*kUu*fkUu/ut;

                scalar uTauNew = ut + f/df;
                err = mag((ut - uTauNew)/ut);
                ut = uTauNew;

            } while (ut > ROOTVSMALL && err > 0.01 && ++iter < 10);

            uTau[faceI] = max(0.0, ut);
        }
    }
    return tuTau;
}

calcUtau 函数，其实是在用牛顿法迭代求解 $y^+$，进而得到 $u_\tau$ 的值。calcNut 函数中
$$
\frac{u_\tau ^2}{|\frac{U_w-U_c}{d}|} - \nu = \frac{\tau_w}{|\frac{U_w-U_c}{d}|} -\nu = \nu_{eff} - \nu = \nu_t
$$

这个壁面函数使用的是从粘性底层连续变化到对数层的 $y^+ \text{-} u^+$ 关系式，所以，这个可以认为是网格无关的，即不管第一层网格落在哪个区，都是有效的。如果网格无法做到全部位于粘性层或者对数区，建议用这个壁面条件。

(5). nutUTabulatedWallFunction

这个壁面函数，需要从外部读取一个 $U^+ \text{-}\,Re_y$ 数据表，通过从这个数据表插值来得到 $U^+$ 的值。其中 $Re_y=yU/\nu$ 。

// 构造函数
nutUTabulatedWallFunctionFvPatchScalarField::
nutUTabulatedWallFunctionFvPatchScalarField
(
    const fvPatch& p,
    const DimensionedField<scalar, volMesh>& iF,
    const dictionary& dict
)
:
    nutWallFunctionFvPatchScalarField(p, iF, dict),
    uPlusTableName_(dict.lookup("uPlusTable")),
    uPlusTable_
    (
        IOobject
        (
            uPlusTableName_,
            patch().boundaryMesh().mesh().time().constant(),
            patch().boundaryMesh().mesh(),
            IOobject::MUST_READ_IF_MODIFIED,
            IOobject::NO_WRITE,
            false
        ),
        true
    )
{}

$U^+$ 和 $\nu_t$ 分别由函数 calcUPlus 和 calcNut 来计算。

tmp<scalarField> nutUTabulatedWallFunctionFvPatchScalarField::calcNut() const
{
    const label patchi = patch().index();
    const turbulenceModel& turbModel =
        db().lookupObject<turbulenceModel>("turbulenceModel");
    const scalarField& y = turbModel.y()[patchi];
    const fvPatchVectorField& Uw = turbModel.U().boundaryField()[patchi];
    const scalarField magUp(mag(Uw.patchInternalField() - Uw));
    const scalarField magGradU(mag(Uw.snGrad()));
    const tmp<volScalarField> tnu = turbModel.nu();
    const volScalarField& nu = tnu();
    const scalarField& nuw = nu.boundaryField()[patchi];
    return
        max
        (
            scalar(0),
            sqr(magUp/(calcUPlus(magUp*y/nuw) + ROOTVSMALL))
           /(magGradU + ROOTVSMALL)
          - nuw
        );
    // magUp/UPlus = utau, sqr(utau) = tauw, tauw/magGradU = nuEff = nut + nu
}

tmp<scalarField> nutUTabulatedWallFunctionFvPatchScalarField::calcUPlus
(
    const scalarField& Rey
) const
{
    tmp<scalarField> tuPlus(new scalarField(patch().size(), 0.0));
    scalarField& uPlus = tuPlus();
    forAll(uPlus, faceI)
    {
        uPlus[faceI] = uPlusTable_.interpolateLog10(Rey[faceI]);
    }
    return tuPlus;
}

注意这里 calcUPlus 用的是 interpolateLog10 函数来插值，这个函数的定义为

template<class Type>
Type Foam::uniformInterpolationTable<Type>::interpolateLog10
(
    scalar x
) const
{
    if (log10_)
    {
        if (x > 0)
        {
            x = ::log10(x);
        }
        else if (bound_ && (x <= 0))
        {
            x = x0_;
        }
        else
        {
            FatalErrorIn
            (
                "uniformInterpolationTable<Type>::interpolateLog10(scalar x)"
            )   << "Table " << name() << nl
                << "Supplied value must be greater than 0 when in log10 mode"
                << nl << "x=" << x << nl << exit(FatalError);
        }
    }
    return interpolate(x); // 这个是普通的线性插值函数
}

即计算 x 的对数（log10），在将计算结果用来进行线性插值。所以，用这个壁面函数的时候，要注意你所提供的数据表是普通线性坐标的还是对数坐标的。

基本上常见的处理壁面上的湍流粘度的方法就是以上几种了。OpenFOAM 中还提供了几个能处理粗糙壁面的壁面函数( nutURoughWallFunction ， nutkRoughWallFunction )，以及处理大气层边界的(nutkAtmRoughWallFunction，需要跟 atmBoundaryLayerInletVelocity 这个入口边界配合使用 )，细节这里不再详述了，有需要时可以去看相关代码，代码结构是类似的，只是具体计算公式不一样。

Giskard's CFD Learning Tricks

CFD and Scientific Computing

OpenFOAM 中的壁面函数（四）

5. 湍流粘度 $\nu_t$ 的壁面函数