[OpenFOAM 소스코드 파해치기] 벽함수 (Wall function)

일본어 원문링크

작성일 : 2012년  9월  15일

번역일 : 2016년  3월   30일

크롬 브라우저로 보시는 것을 권장해 드립니다.

 ------ OpenFOAM 소스코드 파해치기 시리즈 ------

OpenFOAM 소스코드 파해치기 목차로 이동

 ------------------------------------------------------------

 ----- OpenFOAM 소스코드를 다루는 문법 기본 -----

OpenFOAM 소스코드를 다루는 문법 기본으로 이동

 ------------------------------------------------------------

 

들어가기

벽함수 (Wall Function)에 대해 알아보자

사용 버전

OpenFOAM 2.1.1

벽함수

벽 근처에서는 k-epsilon 모델의 가정을 적용할수 없으므로, 근방에서 벽함수라는 특별한 처리를 수행한다. 벽에서부터의 거리를 무차원한 y+, 벽근방의 속도를 무차원한 u+ 를 각각 x축, y축으로 해 그래프로 그리면 선형적으로 증가하는 층류경계인 점성층 (viscous sub-layer, y+ < 11 정도)와 대수적으로 변화하는 대수영역 (30 < y+ < 300 정도)로 구분된다. 

대수영역에서는 u+ 은 다음과 같은 식으로 표현된다.

u+ = (1/κ) ln(E y+)

E, κ 는 매끈한면에서 E = 9.7, κ = 0.42 로 알려져 있다. 이 영역에서는 난류에너지의 생성과 소산이 평행상태이라고 하고, 아래와 같은 가정이 가능하다.

G = ρε

G 는 k 방정식의 생성항으로 다음과 같이 표현된다.

G = μ_t (du/dy)^2

*** 여기서 du/dy는 편미분.

μ_t 는 난류점성계수로 아래과 같이 표현된다.

μ_t = ρ C_μ k^2 / ε

여기서부터 계산하면 아래와 같은 식을 얻을 수 있다.

ε = (C_μ)^(3/4) k^(3/2) / (κy)
G = μt (du/dy) (C_μ)^(1/4) k^(1/2) / (κy)
μ_t = y+/u+

벽면에서 k 의 구배를 0 으로 하여 나머지 값들을 위의 식을 이용해 계산한다.

점성층에서의 u+ 는 아래와 같이 표현된다.

u+ = y+

 

벽함수 경계조건

여기서는 암축성유체의 벽함수 경계조건의 코드를 살펴보자.

$FOAM_SRC/turbulenceModels/compressible/RAS의 derivedFvPatchFields/wallFunction에 있다.

• kqRWallFunctions
• epsilonWallFunctions
• mutWallFunction

kqRWallFunctions는  k 방정식용으로 실제로는 구배 0 의 경계조건을 부여한다.

epsilonWallFunction은 epsilon 방정식용으로 벽근방의 epsilon과 k 방정식의 생성항 G를 계산한다.

epsilonWallFunctionFvPatchScalarField.C

void epsilonWallFunctionFvPatchScalarField::updateCoeffs()
{
    ...

    const RASModel& rasModel = db().lookupObject<RASModel>("RASProperties");

    const scalar Cmu25 = pow025(Cmu_);
    const scalar Cmu75 = pow(Cmu_, 0.75);

    const scalarField& y = rasModel.y()[patchI];

    volScalarField& G =
        const_cast<volScalarField>(db().lookupObject<volScalarField>(GName_));

    DimensionedField<scalar, volMesh>& epsilon =
        const_cast<DimensionedField<scalar, volMesh>&>
        (
            dimensionedInternalField()
        );

    const tmp<volScalarField> tk = rasModel.k();
    const volScalarField& k = tk();

    const scalarField& muw = rasModel.mu().boundaryField()[patchI];

    const tmp<volScalarField> tmut = rasModel.mut();
    const volScalarField& mut = tmut();
    const scalarField& mutw = mut.boundaryField()[patchI];

    const fvPatchVectorField& Uw = rasModel.U().boundaryField()[patchI];

    const scalarField magGradUw(mag(Uw.snGrad()));

    // Set epsilon and G
    forAll(mutw, faceI)
    {
        label faceCellI = patch().faceCells()[faceI];

        epsilon[faceCellI] = Cmu75*pow(k[faceCellI], 1.5)/(kappa_*y[faceI]);

        G[faceCellI] =
            (mutw[faceI] + muw[faceI])
           *magGradUw[faceI]
           *Cmu25*sqrt(k[faceCellI])
           /(kappa_*y[faceI]);
    }

    ...
}

epsilon 은 식에서 보이는 대로, G는 식에서 μ_t 부분이 μ + μ_t 으로 되어있다.

mutWallFunctions는 난류점성의 경계조건을 부여한다.

mutkWallFunctionFvPatchScalarField.C

tmp<scalarField> mutkWallFunctionFvPatchScalarField::calcMut() const
{
    const label patchI = patch().index();
    const RASModel& rasModel = db().lookupObject<RASModel>("RASProperties");
    const scalarField& y = rasModel.y()[patchI];
    const scalarField& rhow = rasModel.rho().boundaryField()[patchI];
    const tmp<volScalarField> tk = rasModel.k();
    const volScalarField& k = tk();
    const scalarField& muw = rasModel.mu().boundaryField()[patchI];

    const scalar Cmu25 = pow025(Cmu_);

    tmp<scalarField> tmutw(new scalarField(patch().size(), 0.0));
    scalarField& mutw = tmutw();

    forAll(mutw, faceI)
    {
        label faceCellI = patch().faceCells()[faceI];

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

        if (yPlus > yPlusLam_)
        {
            mutw[faceI] = muw[faceI]*(yPlus*kappa_/log(E_*yPlus) - 1);
        }
    }

    return tmutw;
}

y+ 가 점성층의 상한 (yPlusLam_) 보다 높은 곳에서  μ_t 를 계산한다. 여기서는 μ_t = y+/u+ 가 아닌 μ + μ_t = y+/u+ 로 계산한다.

yPlusLam_은 다음과 같이 계산된다.

scalar mutkWallFunctionFvPatchScalarField::calcYPlusLam
(
    const scalar kappa,
    const scalar E
) const
{
    scalar ypl = 11.0;

    for (int i=0; i<10; i++)
    {
        ypl = log(E*ypl)/kappa;
    }

    return ypl;
}

점성층의 식 u+ = y+ 를 반복계산으로 풀고 있다.

주의점

이상 설명한 코드에서 μ_t 이외는 대수영역을 전제로 하고 있다. 따라서, 벽면 첫번째 격자의 높이의 y+ 를 대수영역 (30 < y+ < 300 정도)에 넣어야 한다.

추가정보

ANSYS FLUENT에서는 점성층에서 eplison은 대수영역의 식을, G는 점성층의 식을 사용한다.

OpenFOAM 소스코드 파해치기 목차로 이동