[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.);
}

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