본문 바로가기

OpenFOAM/소스코드 파해치기

[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는 점성층의 식을 사용한다.