openfoam學習心得—壁面函數在CFD中如何實施

最近研究壁面函數，發現邱教授部落格，講的很是不錯！記錄一下

http://xiaopingqiu.github.io/2016/04/25/wallFunctions4/

在一般的外流場計算中，一般遇到的邊界有入口、出口、對稱面以及壁面-v7程式設計指南P40對此處應用的邊界條件進行了說明，我們重點關注壁面，壁面一般設定為v-no-slip,p-法向梯度為0，代入離散方程進而進行求解，當加入湍流模型時，還會存在湍流參量k、 ε \varepsilon ε、 ω \omega ω，nut，部落客主要想歸納在運用壁面函數時，到底該如何對這些參量設定邊界條件。

openfoam中的壁面函數按其計算位置分為兩種—第一種計算壁面上的值，第二種計算靠近壁面第一層網格，其中心的值。

openfoam裡面求解器源代碼都位于applications裡面，如壁面函數等預先定義好的類位于src裡面，其主要有六個k、 ε \varepsilon ε、 ω \omega ω，nut、f、v

f、v主要是針對f-v2湍流模型，這裡我們主要看k、 ε \varepsilon ε、 ω \omega ω、nut，注意檔案夾的名稱以s結尾，說明對于某一個參量，可能存在不止一種壁面函數，我們進去下一級目錄檢視，可以發現 ε \varepsilon ε、 ω \omega ω隻有一種壁面函數，k有兩種，nut足足有八種壁面函數可供選擇！今天筆者主要閱讀kqRwallFunctions源代碼，遂做些标注，友善日後查閱也供他人參考，有了解不對的地方，還望大佬不吝賜教！

kqRwallFunction.H-這是類聲明，固沒有出現函數體

namespace Foam         //名稱空間
{
template<class Type>    //模闆類定義符,Type可以是任何類
class kqRWallFunctionFvPatchField  
:
    public zeroGradientFvPatchField<Type> //該類公有繼承于模闆類
{

public:      //公有成員

   
    //- Runtime type information
    TypeName("kqRWallFunction");  //成員變量初始化


    // Constructors

        //- Construct from patch and internal field
        kqRWallFunctionFvPatchField //構造函數 無傳回值，以類名命名，初始化以引用類型傳遞參數
        (
            const fvPatch&,
            const DimensionedField<Type, volMesh>&
        );

        //- Construct from patch, internal field and dictionary
        kqRWallFunctionFvPatchField
        (
            const fvPatch&,
            const DimensionedField<Type, volMesh>&,
            const dictionary&
        ); //構造函數的重載

        //- Construct by mapping given
        // kqRWallFunctionFvPatchField
        //  onto a new patch
        kqRWallFunctionFvPatchField
        (
            const kqRWallFunctionFvPatchField&,
            const fvPatch&,
            const DimensionedField<Type, volMesh>&,
            const fvPatchFieldMapper&
        );

        //- Copy constructor
        kqRWallFunctionFvPatchField
        (
            const kqRWallFunctionFvPatchField&
        );

        //- Construct and return a clone
        virtual tmp<fvPatchField<Type>> clone() const //虛函數，可在子類中重新定義，傳回值是一個tmp<fvPatchField<Type>>類的對象，這是一個模闆類
        {
            return tmp<fvPatchField<Type>>
            (
                new kqRWallFunctionFvPatchField(*this)//*this表示一個kqRWallFunctionFvPatchField 對象，調用了上面的copy構造函數
            );
        }

        //- Copy constructor setting internal field reference
        kqRWallFunctionFvPatchField
        (
            const kqRWallFunctionFvPatchField&,
            const DimensionedField<Type, volMesh>&
        );

        //- Construct and return a clone setting internal field reference
        virtual tmp<fvPatchField<Type>> clone
        (
            const DimensionedField<Type, volMesh>& iF
        ) const
        {
            return tmp<fvPatchField<Type>> //生成tmp<fvPatchField<Type>>類的對象
            (
                new kqRWallFunctionFvPatchField(*this, iF)//其實函數參數是一個指針，指向kqRWallFunctionFvPatchField對象
            );
        }


    // Member Functions---成員函數

        //- Evaluate the patchField
        virtual void evaluate
        (
            const Pstream::commsTypes commsType = Pstream::commsTypes::blocking
        );

        //- Write
        virtual void write(Ostream&) const;
};
           

kqRwallFunction.C

#include "kqRWallFunctionFvPatchField.H"
#include "addToRunTimeSelectionTable.H"

// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //

template<class Type> //同上，模闆類修飾符
Foam::kqRWallFunctionFvPatchField<Type>::kqRWallFunctionFvPatchField
(
    const fvPatch& p,
    const DimensionedField<Type, volMesh>& iF
)
:
    zeroGradientFvPatchField<Type>(p, iF)
{}//調用父類構造函數


template<class Type>
Foam::kqRWallFunctionFvPatchField<Type>::kqRWallFunctionFvPatchField
(
    const fvPatch& p,
    const DimensionedField<Type, volMesh>& iF,
    const dictionary& dict
)
:
    zeroGradientFvPatchField<Type>(p, iF, dict)
{}


template<class Type>
Foam::kqRWallFunctionFvPatchField<Type>::kqRWallFunctionFvPatchField
(
    const kqRWallFunctionFvPatchField& ptf,
    const fvPatch& p,
    const DimensionedField<Type, volMesh>& iF,
    const fvPatchFieldMapper& mapper
)
:
    zeroGradientFvPatchField<Type>(ptf, p, iF, mapper)
{}


template<class Type>
Foam::kqRWallFunctionFvPatchField<Type>::kqRWallFunctionFvPatchField
(
    const kqRWallFunctionFvPatchField& tkqrwfpf
)
:
    zeroGradientFvPatchField<Type>(tkqrwfpf)
{}


template<class Type>
Foam::kqRWallFunctionFvPatchField<Type>::kqRWallFunctionFvPatchField
(
    const kqRWallFunctionFvPatchField& tkqrwfpf,
    const DimensionedField<Type, volMesh>& iF
)
:
    zeroGradientFvPatchField<Type>(tkqrwfpf, iF)
{}
//全是調用父類構造函數，這個類就是個空殼

// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //

template<class Type>
void Foam::kqRWallFunctionFvPatchField<Type>::evaluate
(
    const Pstream::commsTypes commsType
)
{
    zeroGradientFvPatchField<Type>::evaluate(commsType);//調用父類函數
}


template<class Type>
void Foam::kqRWallFunctionFvPatchField<Type>::write(Ostream& os) const
{
    zeroGradientFvPatchField<Type>::write(os); //又是調用父類
    writeEntry(os, "value", *this);
}
           

讀完這份源碼，唯一的感受就是，kqRwallFunction這個邊界條件就是zeroGradient，完全沒有引入新的東西。

kLowReWallFunction.C

這是另外一個k壁面函數，雖然叫做LowRe但是也适用于HighRe，其繼承于fixedValueFvPatchField,可以看出其為一個固定值邊界條件，我們主要分析它的成員函數：

void kLowReWallFunctionFvPatchScalarField::updateCoeffs()
{
    if (updated())//判斷patch單元的值是否更新
    {
        return;
    }

    const label patchi = patch().index();//邊界是被劃分為不同的patch，這是每個patch的标号（label）

    const turbulenceModel& turbModel = db().lookupObject<turbulenceModel>
    (
        IOobject::groupName
        (
            turbulenceModel::propertiesName,
            internalField().group()
        )
    ); //初始化一個對象

    const nutWallFunctionFvPatchScalarField& nutw =
        nutWallFunctionFvPatchScalarField::nutw(turbModel, patchi);

    const scalarField& y = turbModel.y()[patchi]; //算near wall distant
    // tmp<volScalarField>可以了解為volScalarField* 其實也就是個包裝過的指針，搭配new堆記憶體使用
    const tmp<volScalarField> tk = turbModel.k();
    const volScalarField& k = tk();//算第一個單元的k值

    const tmp<scalarField> tnuw = turbModel.nu(patchi);
    const scalarField& nuw = tnuw(); //算第一個單元的粘度值

    const scalar Cmu25 = pow025(nutw.Cmu());//算摩擦速度Cmu25=Cmu^0.25

    scalarField& kw = *this;//初始化修正的k

    // Set k wall values
    forAll(kw, facei)
    {
        label celli = patch().faceCells()[facei];//由面值算單元值

        scalar uTau = Cmu25*sqrt(k[celli]);//算摩擦速度

        scalar yPlus = uTau*y[facei]/nuw[facei];//算y+

        if (yPlus > nutw.yPlusLam())
        {
            scalar Ck = -0.416;
            scalar Bk = 8.366;
            kw[facei] = Ck/nutw.kappa()*log(yPlus) + Bk;//注意這裡是kw[facei]，也就是說邊界上的湍動能是固定值
        }
        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);//個人認為這個時候K場邊界上應該是零梯度的，因為k+應該是面中心的值但是代碼卻是 kw[facei]，說明面中心的值等于面上的值。

    }

    // Limit kw to avoid failure of the turbulence model due to division by kw
    kw = max(kw, small);

    fixedValueFvPatchField<scalar>::updateCoeffs();//執行父類的函數，更新update值
    // TODO: perform averaging for cells sharing more than one boundary face
}