본문 바로가기

OpenFOAM/소스코드 파해치기

[OpenFOAM 소스코드 파해치기] 다공질체 (porous media)

일본어 원문링크


작성일 : 2012년  4월  30일

번역일 : 2016년  3월   29일


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


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

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

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

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

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

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



 

  들어가기 


다공질체 (porous media) 를 다루는 코드를 살펴보자.



  사용 버전


OpenFOAM 2.1.0



  파일  


18-icoFoam.tar.gz 

18-square.tar.gz 



  소스코드 조사  


다공질체를 다루는 솔버로 porousSimpleFoam이 있다. 이것의 UEqn.H 을 살펴보자. 아래와 같은 코드를 통해 다공질체의 저항을 적용한다.

pZones.addResistance(UEqn(), tTU());


pZones는 createPorousZones.H 에서 아래와 같이 정의된다.

porousZones pZones(mesh);


그렇다면, porousZones의 addResistance()를 찾아보자.

$FOAM_SRC/finiteVolume/cfdTools/general/porousMedia/porousZones.H

void Foam::porousZone::addResistance(fvVectorMatrix& UEqn) const { ... const scalarField& V = mesh_.V(); scalarField& Udiag = UEqn.diag(); vectorField& Usource = UEqn.source(); const vectorField& U = UEqn.psi(); .... addViscousInertialResistance ( Udiag, Usource, V, geometricOneField(), mesh_.lookupObject<volScalarField>("nu"), U ); ... }


addViscousInertialResistance() 등으로 실제 저항을 부여하고 있다.

$FOAM_SRC/finiteVolume/cfdTools/general/porousMedia/porousZoneTemplate.H

template<class RhoFIeldType> void Foam::porousZone::addViscousInertialResistance ( scalarField& Udiag, vectorField& Usource, const scalarField& V, const RhoFieldType& rho, const scalarField& mu, const vectorField& U ) const { const tensor& D = D_.value(); const tensor& F = F_.value(); forAll(cellZoneIds_, zoneI) { const labelList& cells = mesh_.cellZones()[cellZoneIds_[zoneI]]; forAll(cells, i) { const tensor dragCoeff = mu[cells[i]]*D + (rho[cells[i]]*mag(U[cells[i]]))*F; const scalar isoDragCoeff = tr(dragCoeff); Udiag[cells[i]] += V[cells[i]]*isoDragCoeff; Usource[cells[i]] -= V[cells[i]]*((dragCoeff - I*isoDragCoeff) & U[cells[i]]); } } }


다공질체의 저항은 보통 운동량방정식의 소스항으로 다루어진다. 위의 코드에서는 U 의 소스 Usource 와 U 의 계수행렬의 대각성분의 값을 설정하고 있다. Usource 에는 저항의 식이 설정되고 있다. Udiag 의 설정은 안정화를 위한 것이다.

 


   소스항에 의한 장애물이 표현 


위에서 살펴본 바와 같이 소스항에 의한 장애물의 표현을 테스트해보자. 운동량방정식에서 속도를 고정하는 단순한 표현은 아래와 같다.

UEqn = A*(U - U0)


A는 큰 값, U0는 고정한 속도이다. A가 충분히 클때, 양변에서 A 를 빼면 U-U0 = 0 즉 U = U0 가 되는 방식이다.


이때, 소스항은 U에 의존하고 있으므로 비선형이며, 불안정할 수 밖에 없다. 따라서, 안정화를 위해 소스항을 다음과 같이 표현한다.

UEqn = dS*U + S0


dS 는 S(U) 의 U 에 의한 미분계수, S0 는 U 의 함수로 볼 수 있다. 이렇게 표현하면, 우변의 첫번째항을 좌변으로 넘길 수 있고, 불안정성이 감소된다. 안정화를 위해 dS는 음수로 하는것이 좋다는 것을 알 수 있다.


위에서 살펴본 것을 바탕으로, 장애물을 다음과 같이 표현할 수 있다.

UEqn = -A*U + A*U0


장애물을 설정하는 코드를 icoFoam에 추가해보자. 먼저, U의 방정식 UEqn에서 계수행렬의 대각성분과 소스항을 취득한다.

scalarField& Udiag = UEqn.diag();

vectorField& Usource = UEqn.source();


A를 설정한다. 여기의 값은 대충 넣어도 된다. 또한, 장애물을 설정하는 영역을 "c0" 라는 셀존으로부터 받도록 한다.

scalar A = 1e8;

sonst labelList& cells = mesh.cellZone()["c0"];


소스항을 설정한다.

forAll(cells, i)

{

Udiag[cells[i]] += A;

}


이때, U0 = 0 이므로, 소스항은 -A*U, 이것을 좌변에 이동시켜 계수행렬에 포함시키므로 위와 같이 표현하게 된다.


영역의 중앙에 셀존 c0을 설정하여, 왼쪽부터 오른쪽으로의 흐름을 만들면 다음과 같이 된다.



U0 에 0 이외의 값을 설정하는 경우, 다음과 같이 할 수 있다.

forAll(cells, i)

{

Udiag[cells[i]] += A;

Usource[cells[i]] += A*vector(0., 1., 0.);

}