2016-04-25

OpenFOAM 中的壁面函数（三）

这篇来看看计算湍动能 $\varepsilon$ 和 $\omega$ 的壁面函数。

3. 湍动能耗散 $\varepsilon$ 的壁面函数

本篇来看看 OpenFOAM 中的 epsilonWallFunction，共有两个： epsilonWallFunction 和 epsilonLowReWallFunction。

(1). epsilonWallFunction

epsilonWallFunction 代码比前面的 kqRWallFunction 复杂多了，主要原因在于这里需要得到的是 epsilon 在临近网格的值，而且，需要考虑包含两个边界面的网格。这里先来梳理代码的脉络，然后再看具体的计算细节。
外部调用的主要是 updateCoeffs() 函数，所以，从这个函数看起。

void epsilonWallFunctionFvPatchScalarField::updateCoeffs()
{
    if (updated())
    {
        return;
    }
    const turbulenceModel& turbulence =
        db().lookupObject<turbulenceModel>(turbulenceModel::typeName);

    setMaster();

    if (patch().index() == master_)
    {
        createAveragingWeights();
        calculateTurbulenceFields(turbulence, G(true), epsilon(true));
    }
    const scalarField& G0 = this->G();
    const scalarField& epsilon0 = this->epsilon();

    typedef DimensionedField<scalar, volMesh> FieldType;

    FieldType& G =
        const_cast<FieldType&>
        (
            db().lookupObject<FieldType>(turbulence.GName())
        );
    //这里是获取内部场，所以，修改这里的引用 "epsilon",相当于修改 epsilon 的内部场值。
    FieldType& epsilon = const_cast<FieldType&>(dimensionedInternalField());

    forAll(*this, faceI)
    {
        label cellI = patch().faceCells()[faceI];

        G[cellI] = G0[cellI];
        epsilon[cellI] = epsilon0[cellI];
    }
    fvPatchField<scalar>::updateCoeffs();
}

一步一步来看。首先是调用了 setMaster() 函数，来看看这个函数以及相关的一个函数 epsilonPatch 的代码：

void epsilonWallFunctionFvPatchScalarField::setMaster()
{
    if (master_ != -1) // 如果当前处理的边界的 master_ != -1，说明它已被处理过，直接返回
    {
        return;
    }
    const volScalarField& epsilon =
        static_cast<const volScalarField&>(this->dimensionedInternalField());

    const volScalarField::GeometricBoundaryField& bf = epsilon.boundaryField();

    label master = -1;
    forAll(bf, patchI)
    {
        if (isA<epsilonWallFunctionFvPatchScalarField>(bf[patchI]))
        {
            epsilonWallFunctionFvPatchScalarField& epf = epsilonPatch(patchI);

            if (master == -1) // 只有头一个被处理的边界满足这个条件
            {
                master = patchI;
            }

            epf.master() = master; // 这意味着所有边界的 master_ 数据成员都将赋值为头一个被处理的边界的编号，即第一个被处理的边界是master
        }
    }
}

epsilonWallFunctionFvPatchScalarField&
epsilonWallFunctionFvPatchScalarField::epsilonPatch(const label patchI)
{
    const volScalarField& epsilon =
        static_cast<const volScalarField&>(this->dimensionedInternalField());

    const volScalarField::GeometricBoundaryField& bf = epsilon.boundaryField();

    const epsilonWallFunctionFvPatchScalarField& epf =
        refCast<const epsilonWallFunctionFvPatchScalarField>(bf[patchI]);

    return const_cast<epsilonWallFunctionFvPatchScalarField&>(epf);
}

从上述代码可以看出， epsilonPatch 函数需要一个参数，这个参数的含义是某一个边界的序号，返回的是指向这个边界的一个 epsilonWallFunctionFvPatchScalarField 类型的引用。
在此基础上，再来看 setMaster。先判断当前边界的数据成员 master_ 是否不等于-1，如果成立则不做任何操作，直接返回；否则，先获取到 epsilon 的所有边界，存在变量 bf 中，然后，遍历 bf ，如果边界的类型是 epsilonWallFunctionFvPatchScalarField，则判断临时变量 master 是否等于 -1，等于则将边界的序号 patchI 赋值给 master，并临时变量 master 的值赋给 patchI 对应边界的数据成员 master_。举个例子，假设有一个算例，有两个边界上使用了 epsilonWallFunctionFvPatchScalarField 类型的边界条件，两个边界的编号分别是 patchI = 0 和 patchI = 1。则在上述循环过程中，当 patchI = 0时， master == -1 肯定成立。于是，patchI = 0 对应边界的数据成员 master_ 被赋值为0；而当遍历到 patchI = 1 时，此时master = 0，所以，结果是 patchI = 1 的边界的数据成员 master_ 也被赋值为0。

继续向下看，如果 patch.index() == master_ ，则调用两个函数。这个怎么理解呢？还以上面的那个简单例子来说明。注意，在外部调用边界条件的时候，也是会依次调用一个场的所有边界的边界条件的。在这里的简单例子中，有两个边界的类型是 epsilonWallFunctionFvPatchScalarField ，所以，我们假设调用 patchI = 0 对应的边界时，由于初始化时数据成员 master_ 赋值为 -1 ，所以，调用 patchI = 0 的边界时， setMaster 函数中的操作会进行。而根据上面的分析，调用 patchI = 0 的边界时， setMaster 函数同时也将 patchI = 1 边界的数据成员 master_ 赋值为 0了，所以，在外部调用 patchI = 1 的边界时， setMaster 函数将不作任何操作，直接返回。同样的，在外部调用 patchI = 0 的边界时，patch.index() == master_ 条件是成立的，所以 createAveragingWeights() 和 calculateTurbulenceFields(turbulence, G(true), epsilon(true)); 两个语句将会执行；而在外部调用 patchI = 1 边界时，由于 patch.index() == master_ 不成立，这两个语句将不执行。

再继续往前看， const scalarField& G0 = this->G(); const scalarField& epsilon0 = this->epsilon(); ，这里是将成员函数 G 和 epsilon 的返回值分别赋给变量 G0 和 epsilon0。开看一下成员函数的定义

scalarField& epsilonWallFunctionFvPatchScalarField::G(bool init)
{
    if (patch().index() == master_) // 只有头一个被处理的边界满足这个条件
    {
        if (init) // init 缺省值是 false 
        {
            G_ = 0.0;
        }
        return G_;
    }
    return epsilonPatch(master_).G(); // 对于不是 master 的边界，返回master边界的数据成员 G_
}

scalarField& epsilonWallFunctionFvPatchScalarField::epsilon(bool init)
{
    if (patch().index() == master_)
    {
        if (init)
        {
            epsilon_ = 0.0;
        }
        return epsilon_;
    }
    return epsilonPatch(master_).epsilon(init);
}

类似的，对于 patchI = 0， patch().index() == master_ ，所以返回值为 patchI = 0 边界的数据成员 G_ 或 epsilon_ (init 的缺省值是 false)；而对于 patchI = 1边界，返回的是 patchI = master_ 对应边界的数据成员 G_ 或 epsilon_，而根据上面的分析， patchI= 1 的边界的数据成员 master_ = 0，因此， patchI = 1 的边界的成员函数返回的是 patchI = 0边界的相应的数据成员。

再往下的内容就很简单了，只是将得到的 G0 和 epsilon0 的值分别赋给当前边界的临近边界网格而已。

到此，代码的框架就基本清晰了，小结一下就是，如果对于某个算例，有多个边界上需要用到 epsilonWallFunctionFvPatchScalarField 类型的边界条件，则，编号更小的那个边界将会被设置成 master。所有的相关计算都在调用 master 边界的时候进行，非 master 的边界，则只需要从 master 那里读取结果即可！

接下来看看外部调用 master 边界的时候，具体做了哪些计算，主要就是看 createAveragingWeights() 和 calculateTurbulenceFields(turbulence, G(true), epsilon(true)); 这两条语句了。

void epsilonWallFunctionFvPatchScalarField::createAveragingWeights()
{
    const volScalarField& epsilon =
        static_cast<const volScalarField&>(this->dimensionedInternalField());

    const volScalarField::GeometricBoundaryField& bf = epsilon.boundaryField();

    const fvMesh& mesh = epsilon.mesh();

    if (initialised_ && !mesh.changing())
    {
        return;
    }

    volScalarField weights
    (
        IOobject
        (
            "weights",
            mesh.time().timeName(),
            mesh,
            IOobject::NO_READ,
            IOobject::NO_WRITE,
            false // do not register
        ),
        mesh,
        dimensionedScalar("zero", dimless, 0.0)
    );

    DynamicList<label> epsilonPatches(bf.size());
    //遍历所有边界，如果边界类型是 epsilonWallFunctionFvPatchScalarField 则将该边界放到 epsilonPatches 这个动态 list 中。
    forAll(bf, patchI)
    {
        if (isA<epsilonWallFunctionFvPatchScalarField>(bf[patchI]))
        {
            epsilonPatches.append(patchI);

            const labelUList& faceCells = bf[patchI].patch().faceCells();
            forAll(faceCells, i)
            {
                label cellI = faceCells[i];
            // weight 衡量的是网格cellI有多少个边界面使用了 epsilonWallFunctionFvPatchScalarField 类型的边界条件
                 weights[cellI]++;
            }
        }
    }
    cornerWeights_.setSize(bf.size());
    
    // 遍历所有 epsilonWallFunctionFvPatchScalarField 类型的边界
    forAll(epsilonPatches, i)
    {
        label patchI = epsilonPatches[i];
        const fvPatchScalarField& wf = weights.boundaryField()[patchI];
    //cornerWeights_存储的所有边界面的weight的倒数，边界面的weight等于其所属网格的weight。所以，如果有一个网格包含两个使用epsilonWallFunction的边界面，那么根据上面的计算，这个网格的weight将是 2，而这两个边界面的 cornerWeights_ 则都是 1/2。 
        cornerWeights_[patchI] = 1.0/wf.patchInternalField();
    }
    // 将数据成员 G_ 和 epsilon_ 初始化为0
    G_.setSize(dimensionedInternalField().size(), 0.0);
    epsilon_.setSize(dimensionedInternalField().size(), 0.0);

    initialised_ = true;
}

void epsilonWallFunctionFvPatchScalarField::calculateTurbulenceFields
(
    const turbulenceModel& turbulence,
    scalarField& G0,
    scalarField& epsilon0
)
{
    // accumulate all of the G and epsilon contributions
    //cornerWeights_ 是一个二维 list，这里是遍历这个list 的第一层
    forAll(cornerWeights_, patchI)
    {
        if (!cornerWeights_[patchI].empty()) // 如果是empty，意味着这个对应的边界不是epsilonWallFunction类型，所以就不需要考虑
        {
            epsilonWallFunctionFvPatchScalarField& epf = epsilonPatch(patchI);

            const List<scalar>& w = cornerWeights_[patchI];

     // 非 empty 则调用 calculate 函数更新 G0 和 epsilon 的值
            epf.calculate(turbulence, w, epf.patch(), G0, epsilon0);
        }
    }
    // apply zero-gradient condition for epsilon
    forAll(cornerWeights_, patchI)
    {
        if (!cornerWeights_[patchI].empty())
        {
            epsilonWallFunctionFvPatchScalarField& epf = epsilonPatch(patchI);

    // 对 epsilon 使用 零梯度边界条件，即将上面计算得到的临近壁面网格的epsilon的值存储在壁面。
            epf == scalarField(epsilon0, epf.patch().faceCells());
        }
    }
}

void epsilonWallFunctionFvPatchScalarField::calculate
(
    const turbulenceModel& turbulence,
    const List<scalar>& cornerWeights,
    const fvPatch& patch,
    scalarField& G,
    scalarField& epsilon
)
{
    const label patchI = patch.index();
    const scalarField& y = turbulence.y()[patchI];
    const scalar Cmu25 = pow025(Cmu_);
    const scalar Cmu75 = pow(Cmu_, 0.75);
    const tmp<volScalarField> tk = turbulence.k();
    const volScalarField& k = tk();
    const tmp<volScalarField> tnu = turbulence.nu();
    const scalarField& nuw = tnu().boundaryField()[patchI];
    const tmp<volScalarField> tnut = turbulence.nut();
    const volScalarField& nut = tnut();
    const scalarField& nutw = nut.boundaryField()[patchI];
    const fvPatchVectorField& Uw = turbulence.U().boundaryField()[patchI];
    const scalarField magGradUw(mag(Uw.snGrad()));

    // Set epsilon and G
    遍历参数 patch 对应的边界的每一个面
    forAll(nutw, faceI)
    {
        label cellI = patch.faceCells()[faceI];
        scalar w = cornerWeights[faceI];
        
        epsilon[cellI] += w*Cmu75*pow(k[cellI], 1.5)/(kappa_*y[faceI]);
        G[cellI] +=
            w
           *(nutw[faceI] + nuw[faceI])
           *magGradUw[faceI]
           *Cmu25*sqrt(k[cellI])
           /(kappa_*y[faceI]);
    }
}

calculate 函数中进行的是实际的计算过程，主要是更新了临近壁面网格的 epsilon 和 G 的值，计算公式如下：
$$
\varepsilon_c = \frac{1}{N} \sum_{f=i}^{N}\left( \frac{c_\mu^{3/4} k_C^{3/2}}{\kappa y_i}\right) \\
\text{相当于} \quad \quad \quad
\varepsilon ^+ = \frac{1}{\kappa y^+} \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad
$$

$$
G_c = \frac{1}{N} \sum_{f=i}^{N}\left( \frac{(\nu + \nu_t)\cdot |\tfrac{U_i-U_c}{d}|\cdot c_\mu^{1/4} k_C^{1/2}}{\kappa y_i}\right)
$$
这里的 Uw.snGrad() 是 fvPatchFields<Type> 类的成员函数：

template<class Type>
Foam::tmp<Foam::Field<Type> > Foam::fvPatchField<Type>::snGrad() const
{
     return patch_.deltaCoeffs()*(*this - patchInternalField());
}

公式中下标 c 表示临近边界的网格， i 表示网格 c 包含的某个边界面元。y 和 d 都表示边界面元所属网格中心到该面元的垂直距离。

还有一个重要的函数， manipulateMatrix

void epsilonWallFunctionFvPatchScalarField::manipulateMatrix
(
    fvMatrix<scalar>& matrix
)
{
    if (manipulatedMatrix())
    {
        return;
    }

    matrix.setValues(patch().faceCells(), patchInternalField());

    fvPatchField<scalar>::manipulateMatrix(matrix);
}

这个函数的功能是修改 matrix 中的值，将当前 patch 每一个面所属网格的值更新到 matrix 中，参考这个帖子。

如果不是使用的低雷诺数湍流模型，则 $\varepsilon$ 应该使用这个边界条件。理论上，边界第一层网格应该设置在对数区。什么是低雷诺数湍流模型呢？这篇帖子的三楼有精彩的解释。

(2). epsilonLowReWallFunction

epsilonLowReWallFunction 继承自 epsilonWallFunction ，在此基础上，增加了一个成员函数 yPlusLam，并重新定义了 calculate 函数

scalar epsilonLowReWallFunctionFvPatchScalarField::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;
}

这个跟 kLowReWallFunction 里是一样的，不再赘述。

void epsilonLowReWallFunctionFvPatchScalarField::calculate
(
    const turbulenceModel& turbulence,
    const List<scalar>& cornerWeights,
    const fvPatch& patch,
    scalarField& G,
    scalarField& epsilon
)
{
    const label patchI = patch.index();
    const scalarField& y = turbulence.y()[patchI];
    const scalar Cmu25 = pow025(Cmu_);
    const scalar Cmu75 = pow(Cmu_, 0.75);
    const tmp<volScalarField> tk = turbulence.k();
    const volScalarField& k = tk();
    const tmp<volScalarField> tnu = turbulence.nu();
    const scalarField& nuw = tnu().boundaryField()[patchI];
    const tmp<volScalarField> tnut = turbulence.nut();
    const volScalarField& nut = tnut();
    const scalarField& nutw = nut.boundaryField()[patchI];
    const fvPatchVectorField& Uw = turbulence.U().boundaryField()[patchI];
    const scalarField magGradUw(mag(Uw.snGrad()));
    
// Set epsilon and G
    forAll(nutw, faceI)
    {
        label cellI = patch.faceCells()[faceI];

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

        scalar w = cornerWeights[faceI];

        if (yPlus > yPlusLam_)
        {
            epsilon[cellI] = w*Cmu75*pow(k[cellI], 1.5)/(kappa_*y[faceI]);
        }
        else
        {
            epsilon[cellI] = w*2.0*k[cellI]*nuw[faceI]/sqr(y[faceI]);
        }
        G[cellI] =
            w
           *(nutw[faceI] + nuw[faceI])
           *magGradUw[faceI]
           *Cmu25*sqrt(k[cellI])
           /(kappa_*y[faceI]);
    }
}

这里需要根据 yPlus 和 yPlusLam_ 的相对大小来选择不同的计算方式。只是，上面这段来自 OpenFOAM-2.3.1 的代码是有问题的！在OpenFOAM-3.0.1 中已经修复成如下

 forAll(nutw, facei)
    {
        label celli = patch.faceCells()[facei];

        scalar yPlus = Cmu25*sqrt(k[celli])*y[facei]/nuw[facei];

        scalar w = cornerWeights[facei];

        if (yPlus > yPlusLam_)
        {
            epsilon0[celli] += w*Cmu75*pow(k[celli], 1.5)/(kappa_*y[facei]);

            G0[celli] +=
                w
               *(nutw[facei] + nuw[facei])
               *magGradUw[facei]
               *Cmu25*sqrt(k[celli])
               /(kappa_*y[facei]);
        }
        else
        {
            epsilon0[celli] += w*2.0*k[celli]*nuw[facei]/sqr(y[facei]);
            G0[celli] += G[celli];
        }
    }
}

yPlus > yPlusLam_ 时，与 epsilonWallFunction 是一样的；
yPlus < yPlusLam_ 时
$$
\varepsilon_c = \frac{1}{N} \sum_{f=i}^{N}\left( \frac{2\cdot k_C \nu_i}{y^2_i}\right)
$$
这个公式等价于
$$
\varepsilon ^+ = 2\frac{k^+}{(y^+)^2}
$$

G 则取在湍流模型中定义的值，不作修改。不过，这里 G0[celli] += G[celli] 意味着假设有一个网格有两个边界面，则这个网格的中计算得到的 G0 ，将是在湍流模型中定义的该网格中的 G 值的 2 倍，即认为每一个边界面对都该网格内的湍动能生成有贡献。

这个边界是给低雷诺数的 $k-\varepsilon$ 模型以及 $v^2\text{-}f$ 模型使用的。用 OpenFOAM-3.0 以下版本的注意了，这些版本的 epsilonLowReWallFunction 有问题，一定不要忘了修正一下上面提到的那个bug ！

4. $\omega$ 的壁面函数

OpenFOAM 中只提供了一个 omegaWallFunction，这个壁面函数，属于一种自动壁面函数，能自动地根据 $y^+$ 的值来在粘性层和对数层切换，过渡层则采用粘性层和对数层混合的结果。
omegaWallFunction 与 epsilonWallFunction 类似，也是需要计算 $\omega$ 和 $P_k$ 在临近边界网格里的值，因此也需要考虑一个网格包含两个以上边界面的情况。具体处理方法跟 epsilonWallFunction 是一样的，所以这里就不重复了，只看具体的计算 $\omega$ 和 $P_k$ 的公式

void omegaWallFunctionFvPatchScalarField::calculate
(
    const turbulenceModel& turbulence,
    const List<scalar>& cornerWeights,
    const fvPatch& patch,
    scalarField& G,
    scalarField& omega
)
{
    const label patchI = patch.index();
    const scalarField& y = turbulence.y()[patchI];
    const scalar Cmu25 = pow025(Cmu_);
    const tmp<volScalarField> tk = turbulence.k();
    const volScalarField& k = tk();
    const tmp<volScalarField> tnu = turbulence.nu();
    const scalarField& nuw = tnu().boundaryField()[patchI];
    const tmp<volScalarField> tnut = turbulence.nut();
    const volScalarField& nut = tnut();
    const scalarField& nutw = nut.boundaryField()[patchI];
    const fvPatchVectorField& Uw = turbulence.U().boundaryField()[patchI];
    const scalarField magGradUw(mag(Uw.snGrad()));

    // Set omega and G
    forAll(nutw, faceI)
    {
        label cellI = patch.faceCells()[faceI];
        scalar w = cornerWeights[faceI];
        scalar omegaVis = 6.0*nuw[faceI]/(beta1_*sqr(y[faceI]));
        scalar omegaLog = sqrt(k[cellI])/(Cmu25*kappa_*y[faceI]);
        omega[cellI] += w*sqrt(sqr(omegaVis) + sqr(omegaLog));
        G[cellI] +=
            w
           *(nutw[faceI] + nuw[faceI])
           *magGradUw[faceI]
           *Cmu25*sqrt(k[cellI])
           /(kappa_*y[faceI]);
    }
}

这里， omegaVis 和 omegaLog 分别指的是在假定第一层网格位于粘性底层和对数层时得到的 omega 的解析解
$$
\omega_{Vis} = \frac{6.0\nu}{\beta_1y^2} \
\omega_{Log} = \frac{k_C^{1/2}}{C_\mu^{1/4}\kappa y}
$$
然后，将 $\omega_{Vis}$ 和 $\omega_{Log}$ 用一个函数混合起来，就得到了
$$
\omega = \sqrt{\omega_{Vis}^2 + \omega_{Log}^2}
$$
只是，这里的湍动能生成项，却似乎并没有使用混合的方法，而是用的基于对数律的公式：
$$
G = \frac{(\nu + \nu_t)\cdot |\frac{U_c-U_w}{d}|\cdot C_\mu^{1/4}k_C^{1/2}}{\kappa y}
$$

$omega$ 方程是能直接积分到壁面，所以，如果使用基于 $\omega$ 的湍流模型，$\omega$ 变量直接使用这个边界条件就可以了。

Giskard's CFD Learning Tricks

CFD and Scientific Computing

OpenFOAM 中的壁面函数（三）

3. 湍动能耗散 $\varepsilon$ 的壁面函数

4. $\omega$ 的壁面函数