OpenFOAM 中的壁面函数(四)

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

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

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

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

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

1
2
3
4
5
6
7
8
9
10
11
void nutWallFunctionFvPatchScalarField::updateCoeffs()
{
if (updated())
{
return;
}

operator==(calcNut());

fixedValueFvPatchScalarField::updateCoeffs();
}

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

  • (1). nutkWallFunction
    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    15
    16
    17
    18
    19
    20
    21
    22
    23
    24
    25
    26
    27
    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^+$。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
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
    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    15
    16
    17
    18
    19
    20
    21
    22
    23
    24
    25
    26
    27
    28
    29
    30
    31
    32
    33
    34
    35
    36
    37
    38
    39
    40
    41
    42
    43
    44
    45
    46
    47
    48
    49
    50
    51
    52
    53
    54
    55
    56
    57
    58
    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。
    1
    2
    3
    4
    tmp<scalarField> nutLowReWallFunctionFvPatchScalarField::calcNut() const
    {
    return tmp<scalarField>(new scalarField(patch().size(), 0.0))
    ;

    }

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

1
2
3
4
5
6
7
8
9
10
11
12
13
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]
    $$
    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    15
    16
    17
    18
    19
    20
    21
    22
    23
    24
    25
    26
    27
    28
    29
    30
    31
    32
    33
    34
    35
    36
    37
    38
    39
    40
    41
    42
    43
    44
    45
    46
    47
    48
    49
    50
    51
    52
    53
    54
    55
    56
    57
    58
    59
    60
    61
    62
    63
    64
    65
    66
    67
    68
    69
    70
    71
    72
    73
    74
    75
    76
    77
    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$ 。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
// 构造函数
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$ 分别由函数 calcUPluscalcNut 来计算。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
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 函数来插值,这个函数的定义为

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
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 中还提供了几个能处理粗糙壁面的壁面函数( nutURoughWallFunctionnutkRoughWallFunction ),以及处理大气层边界的(nutkAtmRoughWallFunction,需要跟 atmBoundaryLayerInletVelocity 这个入口边界配合使用 ),细节这里不再详述了,有需要时可以去看相关代码,代码结构是类似的,只是具体计算公式不一样。