[OpenFOAM 소스코드 파해치기] interFoam

일본어 원문링크

작성일 : 2012년 11월  11일

번역일 : 2016년  4월   5일

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

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

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

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

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

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

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

 

들어가기

interFoam 의 코드를 살펴보자

사용 버전

OpenFOAM 2.1.1

Volume of Fluid (VOF) 법

interFoam은 Volume Of Fluid (VOF) 법에 의한 다상유체(2개의 서로 다른 상) 솔버이다. VOF 법 에서는, 상의 구별에 각각의 상의 체적분율 (0~1) 을 이용한다. 2개의 상이일때 첫번째의 상의 체적분율을 α1 이라고 하면, 두번째 상의 체적분율은 1- α1 으로 계산하게 된다.

상의 체적분율의 거동은, 수송방정식으로 표현한다.

Dα1/Dt = 0

운동량방정식은 단상일때의 식과 같은 것을 2개 풀지만 밀도와 점성을 두개의 유체 각각의 것으로부터 구해낸다.

ρ = α1 ρ1 + (1 - α1) ρ2
μ = α1 μ1 + (1 - α1) μ2

운동량방적식에는 표면장력을 추가한다.

fs = σ κ n  (표면장력)
n = ∇α1/|∇α1|  (경계면의 법선방향 벡터)
κ = ∇・n  (경계면의 곡률)

 

Interface Compression

interFoam에서는, 상의 체적분율을 수송방정식으로 다음과 같이 풀어낸다.

Dα1/Dt + ∇・(α1 (1 - α1) Ur) = 0

Ur 은 상의 상대속도 (U1-U2)로, "compression velocity"라고 불린다. Ur 에 대응하는 면의 유체 플럭스는 다음과 같이 계산된다.

Φr = nf min(Ca |Φ|/|S|, max(|Φ|/|S|))

Ca 는 압축 (compression) 의 정도를 나타내는 정수이다. 이것은 fvSolution에서 지정ㅇ하는 cAlpha에 해당한다.

compression velocity를 이산화하기위해 "interfaceCompression" 이라는 스킴이 사용된다.

fvSchemes

divSchemes
{
...
div(phi,alpha) Gauss vanLeer;
div(phirb,alpha) Gauss interfaceCompression;
}

 

Interface Compression

코드를 살펴보자.

interFoam.C

    while (runTime.run())
    {
        #include "readTimeControls.H"
        #include "CourantNo.H"
        #include "alphaCourantNo.H"
        #include "setDeltaT.H"

        runTime++;

        Info<< "Time = " << runTime.timeName() << nl << endl;

        twoPhaseProperties.correct();

        #include "alphaEqnSubCycle.H"

        // --- Pressure-velocity PIMPLE corrector loop
        while (pimple.loop())
        {
            #include "UEqn.H"

            // --- Pressure corrector loop
            while (pimple.correct())
            {
                #include "pEqn.H"
            }

            if (pimple.turbCorr())
            {
                turbulence->correct();
            }
        }

        runTime.write();

        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
            << nl << endl;
    }

twoPhaseMixture

먼저, twoPhaseProperties.correct() 가 호출되고 있다. twoPhaseProperties 는 createFields.H 에 다음과 같이 정의된다.

twoPhaseMixture twoPhaseProperties(U, phi);

twoPhaseMixture의 correct() 는 다음과 같이 되어있다.

$FOAM_SRC/transportModels/incompressible/incompressibleTwoPhaseMixture/twoPhaseMixture.H

virtual void correct()
{
calcNu();
}

calcNu() 는 다음과 같다. (twoPhaseMixture.C).

//- Calculate and return the laminar viscosity
void Foam::twoPhaseMixture::calcNu()
{
    nuModel1_->correct();
    nuModel2_->correct();

    const volScalarField limitedAlpha1
    (
        "limitedAlpha1",
        min(max(alpha1_, scalar(0)), scalar(1))
    );

    // Average kinematic viscosity calculated from dynamic viscosity
    nu_ = mu()/(limitedAlpha1*rho1_ + (scalar(1) - limitedAlpha1)*rho2_);
}

동점성계수를 계산하고 있다.

alphaEqnSubCycle

다음으로, 다음의 행을 살펴보자.

#include "alphaEqnSubCycle.H"

alphaEqnSubCycle.H

label nAlphaCorr(readLabel(pimple.dict().lookup("nAlphaCorr")));

label nAlphaSubCycles(readLabel(pimple.dict().lookup("nAlphaSubCycles")));

if (nAlphaSubCycles > 1)
{
    dimensionedScalar totalDeltaT = runTime.deltaT();
    surfaceScalarField rhoPhiSum(0.0*rhoPhi);

    for
    (
        subCycle<volScalarField> alphaSubCycle(alpha1, nAlphaSubCycles);
        !(++alphaSubCycle).end();
    )
    {
        #include "alphaEqn.H"
        rhoPhiSum += (runTime.deltaT()/totalDeltaT)*rhoPhi;
    }

    rhoPhi = rhoPhiSum;
}
else
{
    #include "alphaEqn.H"
}

interface.correct();

rho == alpha1*rho1 + (scalar(1) - alpha1)*rho2;

nAlphaSubCycle (fvSolution 에서 지정) 이 1 이상이면, 서브사이클에 들어가 "alphaEqn.H" 로 들어간다. 아니면 서브사이클 없이 진행한다.

subCycle

subCycle을 쫒아가 보자.

$FOAM_SRC/OpenFOAM/algorithms/subCycle/subCycle.H

        subCycle(GeometricField& gf, const label nSubCycles)
        :
            subCycleField<GeometricField>(gf),
            subCycleTime(const_cast<Time&>(gf.time()), nSubCycles)
        {}

subCycleField 는 받은 필드를 유지한다.

$FOAM_SRC/OpenFOAM/db/Time/subCycleTime.C

Foam::subCycleTime::subCycleTime(Time& t, const label nSubCycles)
:
    time_(t),
    nSubCycles_(nSubCycles),
    subCycleIndex_(0)
{
    time_.subCycle(nSubCycles_);
}

최종적으로 Time 까지 수행하고 있다.

Time.C

Foam::TimeState Foam::Time::subCycle(const label nSubCycles)
{
    subCycling_ = true;
    prevTimeState_.set(new TimeState(*this));

    setTime(*this - deltaT(), (timeIndex() - 1)*nSubCycles);
    deltaT_ /= nSubCycles;
    deltaT0_ /= nSubCycles;
    deltaTSave_ = deltaT0_;

    return prevTimeState();
}


void Foam::Time::endSubCycle()
{
    if (subCycling_)
    {
        subCycling_ = false;
        TimeState::operator=(prevTimeState());
        prevTimeState_.clear();
    }
}

서브사이클용으로 시각을 더욱 잘게 나누거나 하는 작업이 수행된다. 서브사이클이 끝나면 원래로 돌려놓는다. 즉, 서브사이클 내에서 보통 시간을 다룰때 서브사이클내부시간으로 다루어진다는 것을 알 수 있다.

alphaEqn

자, alphaEqn.H 으로 가보자.

alphaEqn.H

{
    word alphaScheme("div(phi,alpha)");
    word alpharScheme("div(phirb,alpha)");

    surfaceScalarField phic(mag(phi/mesh.magSf()));
    phic = min(interface.cAlpha()*phic, max(phic));
    surfaceScalarField phir(phic*interface.nHatf());

    for (int aCorr=0; aCorr<nAlphaCorr; aCorr++)
{
         surfaceScalarField phiAlpha
        (
            fvc::flux
            (
                phi,
                alpha1,
                alphaScheme
            )
          + fvc::flux
            (
                -fvc::flux(-phir, scalar(1) - alpha1, alpharScheme),
                alpha1,
                alpharScheme
            )
        );

        MULES::explicitSolve(alpha1, phi, phiAlpha, 1, 0);

        rhoPhi = phiAlpha*(rho1 - rho2) + phi*rho2;
}

Info<< "Phase-1 volume fraction = "
    << alpha1.weightedAverage(mesh.Vsc()).value()
    << "  Min(alpha1) = " << min(alpha1).value()
    << "  Max(alpha1) = " << max(alpha1).value()     
 << endl;
}

여기서는 alpha1 의 수송방정식을 풀고 있다. aphir 이 compression velocity 의 유체플럭스이다. 루프에 들어가서 방정식을 풀고 있는데, nAlphaCorr > 1 이어야 한다 (fvSolution에서 지정). phiAlpha를 만들어 MULES::explicitSolve()로 풀고 있다. phiAlpha의 첫번째 항이 유체이동항, 두번째 항이 compression velocity의 항이다. 두번째항에서 alpharScheme("div(phirb,alpha)") 가 사용되고 있따. 이것은 "interfaceCompression" 지정할 필요가 있음을 의미한다 (fvSchemes).

interfaceCompression

interfaceCompression 스킴은, 다음과 같이 되어 있다.

$FOAM_SRC/transportModels/interfaceProperties/interfaceCompression/interfaceCompression.H

class interfaceCompressionLimiter
{

public:

    interfaceCompressionLimiter(Istream&)
    {}

    scalar limiter
    (
        const scalar cdWeight,
        const scalar faceFlux,
        const scalar phiP,
        const scalar phiN,
        const vector&,
        const scalar
    ) const
    {
        // Quadratic compression scheme
        //return min(max(4*min(phiP*(1 - phiP), phiN*(1 - phiN)), 0), 1);

        // Quartic compression scheme
        return
            min(max(
            1 - max(sqr(1 - 4*phiP*(1 - phiP)), sqr(1 - 4*phiN*(1 - phiN))),
            0), 1);
    }
};

잘 이해가 안되므로, 일단 다음으로 넘어가 보자.

MULES

MULES 는, MULES.H 에 아래와 같이 설명이 있다.

MULES: Multidimensional universal limiter with explicit solution.

Solve a convective-only transport equation using an explicit universal
multi-dimensional limiter.

Parameters are the variable to solve, the normal convective flux and the
actual explicit flux of the variable which is also used to return limited
flux used in the bounded-solution.

.... 다음으로 넘어가자.

$FOAM_SRC/finiteVolume/fvMatrices/solvers/MULES/MULES.H

    surfaceScalarField muEff
    (
        "muEff",
        twoPhaseProperties.muf()
      + fvc::interpolate(rho*turbulence->nut())
    );

    fvVectorMatrix UEqn
    (
        fvm::ddt(rho, U)
      + fvm::div(rhoPhi, U)
      - fvm::laplacian(muEff, U)
      - (fvc::grad(U) & fvc::grad(muEff))
    //- fvc::div(muEff*(fvc::interpolate(dev(fvc::grad(U))) & mesh.Sf()))
    );

    UEqn.relax();

    if (pimple.momentumPredictor())
    {
        solve
        (
            UEqn
         ==
            fvc::reconstruct
            (
                (
                    fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)
                  - ghf*fvc::snGrad(rho)
                  - fvc::snGrad(p_rgh)
                ) * mesh.magSf()
            )
        );
    }

twoPhaseProperties.muf() 는 아래와 같이 계산된다.

twoPhaseMixture.C

Foam::tmp<Foam::surfaceScalarField> Foam::twoPhaseMixture::muf() const
{
    const surfaceScalarField alpha1f
    (
        min(max(fvc::interpolate(alpha1_), scalar(0)), scalar(1))
    );

    return tmp
    (
        new surfaceScalarField
        (
            "muf",
            alpha1f*rho1_*fvc::interpolate(nuModel1_->nu())
          + (scalar(1) - alpha1f)*rho2_*fvc::interpolate(nuModel2_->nu())
        )
    );
}

표면장력

운동량방정식의 우변이 복잡하게 표현되어 있는데, 표면장력이 추가되었기 때문이다. 아래의 항이 표면장력이다.

fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)

interface 는 createFields.H 로 정의되고 있다.

// Construct interface from alpha1 distribution
interfaceProperties interface(alpha1, U, twoPhaseProperties);

$FOAM_SRC/transportModels/interfaceProperties/interfaceProerties.H

tmp<volScalarField> sigmaK() const
{
return sigma_*K_;
}

sigmaK() 는 σκ 를 반납한다. κ 는 아래와 같이 계산되고 있다.

interfaceProperties.C

void Foam::interfaceProperties::calculateK()
{
    const fvMesh& mesh = alpha1_.mesh();
    const surfaceVectorField& Sf = mesh.Sf();

    // Cell gradient of alpha
    const volVectorField gradAlpha(fvc::grad(alpha1_));

    // Interpolated face-gradient of alpha
    surfaceVectorField gradAlphaf(fvc::interpolate(gradAlpha));

    //gradAlphaf -=
    //    (mesh.Sf()/mesh.magSf())
    //   *(fvc::snGrad(alpha1_) - (mesh.Sf() & gradAlphaf)/mesh.magSf());

    // Face unit interface normal
    surfaceVectorField nHatfv(gradAlphaf/(mag(gradAlphaf) + deltaN_));
    correctContactAngle(nHatfv.boundaryField(), gradAlphaf.boundaryField());

    // Face unit interface normal flux
    nHatf_ = nHatfv & Sf;

    // Simple expression for curvature
    K_ = -fvc::div(nHatf_);

    // Complex expression for curvature.
    // Correction is formally zero but numerically non-zero.
    /*
    volVectorField nHat(gradAlpha/(mag(gradAlpha) + deltaN_));
    forAll(nHat.boundaryField(), patchi)
    {
        nHat.boundaryField()[patchi] = nHatfv.boundaryField()[patchi];
    }

    K_ = -fvc::div(nHatf_) + (nHat & fvc::grad(nHatfv) & nHat);
    */
}

중력

중력의 항은 아래와 같다.

- ghf*fvc::snGrad(rho)

ghf 는 createFields.H 에 다음과 같이 정의도어 있다.

surfaceScalarField ghf("ghf", g & mesh.Cf());

중력가속도벡터와 면의 중심위치 벡터의 내적을 나타내고 있다. 즉 중력방향 높이와 중력가속도의 곱을 계산하고 있다.

pEqn.H

pEqn.H 로 가보자.

pEqn.H

{
    volScalarField rAU(1.0/UEqn.A());
    surfaceScalarField rAUf(fvc::interpolate(rAU));

    U = rAU*UEqn.H();
    surfaceScalarField phiU
    (
        "phiU",
        (fvc::interpolate(U) & mesh.Sf())
      + fvc::ddtPhiCorr(rAU, rho, U, phi)
    );

    adjustPhi(phiU, U, p_rgh);

    phi = phiU +
    (
        fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)
      - ghf*fvc::snGrad(rho)
    )*rAUf*mesh.magSf();

    while (pimple.correctNonOrthogonal())
    {
        fvScalarMatrix p_rghEqn
        (
            fvm::laplacian(rAUf, p_rgh) == fvc::div(phi)
        );

        p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell));

        p_rghEqn.solve(mesh.solver(p_rgh.select(pimple.finalInnerIter())));

        if (pimple.finalNonOrthogonalIter())
        {
            phi -= p_rghEqn.flux();
        }
    }

    U += rAU*fvc::reconstruct((phi - phiU)/rAUf);
    U.correctBoundaryConditions();

    #include "continuityErrs.H"

    p == p_rgh + rho*gh;

    if (p_rgh.needReference())
    {
        p += dimensionedScalar
        (
            "p",
            p.dimensions(),
            pRefValue - getRefCellValue(p, pRefCell)
        );
        p_rgh = p - rho*gh;
    }
}

여기서도 표면장력을 고려하고 있다.

압력은 p_rgh 를 풀고 있으며, 이것은 다음과 같이 계산된다.

p_rgh = p - rgh*gh;

gh 는 createFields.H 에 있다.

volScalarField gh("gh", g * mesh.C());

ghf 와 같다. 여기서는 셀의 중심을 사용한다. 즉 p_rgh 는 p- ρgh 를 의미한다.

 

참고문헌

H. Rusche, Computational fluid dynamics of dispersed two-phase flows at high phase fractions, PhD Thesis, Imperial College, 2002

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