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

일본어 원문링크

작성일 : 2011년  12월    3일
번역일 : 2016년   3월  28일

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

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

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

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

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

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

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

 

들어가기

simpleFoam의 코드를 분석해 보자.

사용 버전

OpenFOAM 2.0.1

 

simpleFoam의 파일구성

$FOAM_SOLVERS/incompressible/simpleFoam

• Make/files
• Make/options
• simpleFoam.C
• createFields.H
• UEqn.H
• pEqn.H
• 그외 파생 솔버

simpleFoam의 개요

 

simpleFoam은 SIMPLE법을 이용한 비압축성유체의 비정상난류해석솔버이다. 유체의 방정식은 우동방정식과 연속방정식으로 구성되며, SIMPLE법은 그것을 운동량방정식과 압력수정방정식의 2개의 방정식의 형태로 푼다. 한번에 풀수 없으므로(연속방정식이 만족되지 않는다) 여러번의 반복계산을 통해 해를 얻어낸다. 수렴조건에 만적하면(방정식의 Residual이 작아지면) 계산을 종료하게 된다.

 

SIMPLE법

정상상태의 비압축성유체의 방정식은 아래와 같다.

div(phi u) = -grad(p) + (viscous term) ... (1)

div(u) = 0 ... (2)

식 (1)은 운동량방정식, 식 (2)는 연속방정식이다. 운동량방정식의 부분을 이산화 하여 다음과 같이 나타낸다.

A u = H - grad(p) ... (3)

식 (3)을 u 에 대한 식으로 고치면,

u = H/A - 1/A * grad(p) ... (4)

식 (4)에 divergence를 취해 연속방정식을 만든다.

div(u) = div(H/A) - 1/A * laplacian(p) = 0 ... (5)

식 (5)에 의해, 아래와 같은 압력수정방정식이 만들어진다.

1/A * laplacian(p) = div(H/A) ... (6)

자세한 설명은 생략하도록 하고, SIMPLE법은 아래와 같은 순서로 계산을 수행한다.

1. 운동량방정식 (식 (1))을 푼다.

2. 압력수정방정식 (식 (6))을 푼다.

3. 얻어진 압력으로 식 (4) 에 의해 속도를 수정한다.

4. 위 3과정을 반복한다.

 

simpleFoam의 코드

simpleFoam의 코드를 살펴보자.

simpleFoam.C

#include "fvCFD.H"
#include "singlePhaseTransportModel.H"
#include "RASModel.H"
#include "simpleControl.H"

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

int main(int argc, char *argv[])
{
    #include "setRootCase.H"
    #include "createTime.H"
    #include "createMesh.H"
    #include "createFields.H"
    #include "initContinuityErrs.H"

    simpleControl simple(mesh);

    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

    Info<< "\nStarting time loop\n" << endl;

    while (simple.loop())
    {
        Info<< "Time = " << runTime.timeName() << nl << endl;

        p.storePrevIter();

        // --- Pressure-velocity SIMPLE corrector
        {
            #include "UEqn.H"
            #include "pEqn.H"
        }

        turbulence->correct();

        runTime.write();

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

    Info<< "End\n" << endl;

    return 0;
}

코드초반의 해석환경 설정 후 SIMPLE법의 루프에 들어가서 운동량방정식과 압력수정방정식을 풀어내고 있다. 해석환경 설정의 createFields.H 부터 살펴보도록 하자.

createFields.H

    Info<< "Reading field p\n" << endl;
    volScalarField p
    (
        IOobject
        (
            "p",
            runTime.timeName(),
            mesh,
            IOobject::MUST_READ,
            IOobject::AUTO_WRITE
        ),
        mesh
    );

    Info<< "Reading field U\n" << endl;
    volVectorField U
    (
        IOobject
        (
            "U",
            runTime.timeName(),
            mesh,
            IOobject::MUST_READ,
            IOobject::AUTO_WRITE
        ),
        mesh
    );

    #include "createPhi.H"


    label pRefCell = 0;
    scalar pRefValue = 0.0;
    setRefCell(p, mesh.solutionDict().subDict("SIMPLE"), pRefCell, pRefValue);

    singlePhaseTransportModel laminarTransport(U, phi);

    autoPtr turbulence
    (
        incompressible::RASModel::New(U, phi, laminarTransport)
    );

처음부분에 속도장 U와 압력장 p를 준비시킨다. 이후 createPhi.H 는 아래와 같이 되어있다.

$FOAM_SRC/finiteVolume/cfdTools/incompressible/createPhi.H

surfaceScalarField phi
(
    IOobject
    (
        "phi",
        runTime.timeName(),
        mesh,
        IOobject::READ_IF_PRESENT,
        IOobject::AUTO_WRITE
    ),
    linearInterpolate(U) & mesh.Sf()
);

phi는 페이스(면)의 속도와 면적벡터의 내적이다.

    singlePhaseTransportModel laminarTransport(U, phi);

여기서 constant/transportProperties 를 불러오고 있다.

    autoPtr turbulence
    (
        incompressible::RASModel::New(U, phi, laminarTransport)
    );

여기서 constant/RASProperties 를 불러온다.

simpleFoam.C 로 돌아가보자.

    simpleControl simple(mesh);

이것은, 이전의 버전에서는 system/fvSolution에서부터 하나씩 불러오던  SIMPLE의 설정을 한번에 불러오도록 하는 부분에 해당한다. 수렴판정도 여기서 불러오게 된다.

    while (simple.loop())
{
...
}

SIMPLE법의 반복계산 루프이다.

    p.storePrevIter();

GeometricField::storePrevIter()는, 현재의 값을 나중에 사용하도록 저장한다. 여기서는 이후에 나오는 "p.relax()"를 위해 사용되고 있다.

   {
#include "UEqn.H"
#include "pEqn.H"
}

운동량방정식과 압력수정방정식을 풀고 있다.

UEqn.H

    tmp<fvVectorMatrix> UEqn
    (
        fvm::div(phi, U)
      + turbulence->divDevReff(U)
    );

    UEqn().relax();

    solve(UEqn() == -fvc::grad(p));

운동량방정식을 만들고 풀어내는 부분이다. UEqn은 운동량방정식에서 압력구배항을 제외한 식이다. 이것은 위에서 기술한 SIMPLE법의 순서 중에 "H/A(momentum predictor)"를 계산할 필요가 있기 때문이다.

turbulence->divDevReff(U)는 점성응력의 div취한 항(편차응력항)이다. 다음과 같이 계산된다. 여기서는 k-epsilon 모델의 것을 살펴보도록 하자.

$FOAM_SRC/turbulneceModels/incompressible/RAS/kEpsilon/kEpsilon.C

tmp<fvVectorMatrix> kEpsilon::divDevReff(volVectorField& U) const
{
    return
    (
      - fvm::laplacian(nuEff(), U)
      - fvc::div(nuEff()*dev(T(fvc::grad(U))))
    );
}

속도구배의 전치 항은 양함수 형태로 다루어지고 있다.

UEqn().relax();

relaxation coefficient를 방정식의 계수행렬에 적용시키고 있다.

solve(UEqn() == -fvc::grad(p));

운동량방정식을 풀고 있다. 압력구방행은 양함수 형태로 다루어진다.

압력수정방정식을 살펴보자.

pEqn.C

    p.boundaryField().updateCoeffs();

    volScalarField rAU(1.0/UEqn().A());
    U = rAU*UEqn().H();
    UEqn.clear();

    phi = fvc::interpolate(U, "interpolate(HbyA)") & mesh.Sf();
    adjustPhi(phi, U, p);

    // Non-orthogonal pressure corrector loop
    for (int nonOrth=0; nonOrth<=simple.nNonOrthCorr(); nonOrth++)
    {
        fvScalarMatrix pEqn
        (
            fvm::laplacian(rAU, p) == fvc::div(phi)
        );

        pEqn.setReference(pRefCell, pRefValue);

        pEqn.solve();

        if (nonOrth == simple.nNonOrthCorr())
        {
            phi -= pEqn.flux();
        }
    }

    #include "continuityErrs.H"

    // Explicitly relax pressure for momentum corrector
    p.relax();

    // Momentum corrector
    U -= rAU*fvc::grad(p);
    U.correctBoundaryConditions();

압력수정방정식을 구성해 풀고 있다. 격자의 뒤틀림을 보정하기 위한 반복계산을 포함한다.

p.boundaryFields().updateCoeffs();

p의 경계조건에 대해 계산하고 있다.

    volScalarField rAU(1.0/UEqn().A());
    U = rAU*UEqn().H();
    UEqn.clear();

U = H/A (momentum predictor)를 계산하고 있다.

    phi = fvc::interpolate(U, "interpolate(HbyA)") & mesh.Sf();
    adjustPhi(phi, U, p);

phi 를 계산한다. fvc::interpolate() 로 U의 보간값을 계산하여, 면적벡터와 내적한다.

$FOAM_SRC/finiteVolume/cfdTools/general/adjustPhi/adjustPhi.C

bool Foam::adjustPhi
(
    surfaceScalarField& phi,
    const volVectorField& U,
    volScalarField& p
)
{
...

        scalar massIn = 0.0;
        scalar fixedMassOut = 0.0;
        scalar adjustableMassOut = 0.0;

        forAll(phi.boundaryField(), patchi)
        {
            const fvPatchVectorField& Up = U.boundaryField()[patchi];
            const fvsPatchScalarField& phip = phi.boundaryField()[patchi];

            if (!isA<processorFvsPatchScalarField>(phip))
            {
                if
                (
                    Up.fixesValue()
                 && !isA<inletOutletFvPatchVectorField>(Up)
                )
                {
                    forAll(phip, i)
                    {
                        if (phip[i] < 0.0)
                        {
                            massIn -= phip[i];
                        }
                        else
                        {
                            fixedMassOut += phip[i];
                        }
                    }
                }
                else
                {
                    forAll(phip, i)
                    {
                        if (phip[i] < 0.0)
                        {
                            massIn -= phip[i];
                        }
                        else
                        {
                            adjustableMassOut += phip[i];
                        }
                    }
                }
            }
        }

        ...

        scalar massCorr = 1.0;
        scalar magAdjustableMassOut = mag(adjustableMassOut);

        ...
	
            massCorr = (massIn - fixedMassOut)/adjustableMassOut;
        
	...

        forAll(phi.boundaryField(), patchi)
        {
            const fvPatchVectorField& Up = U.boundaryField()[patchi];
            fvsPatchScalarField& phip = phi.boundaryField()[patchi];

            if (!isA<processorFvsPatchScaclarField>(phip))
            {
                if
                (
                    !Up.fixesValue()
                 || isA<inletOutletFvPatchVectorField>(Up)
                )
                {
                    forAll(phip, i)
                    {
                        if (phip[i] > 0.0)
                        {
                            phip[i] *= massCorr;
                        }
                    }
                }
            }
        }

	...
}

adjustPhi()는 유속에 의해 결정되는 경계의 유체 flux를 질량이 보전되도록 조정한다.

    for (int nonOrth=0; nonOrth<=simple.nNonOrthCorr(); nonOrth++)
    {
        fvScalarMatrix pEqn
        (
            fvm::laplacian(rAU, p) == fvc::div(phi)
        );

        pEqn.setReference(pRefCell, pRefValue);

        pEqn.solve();

        if (nonOrth == simple.nNonOrthCorr())
        {
            phi -= pEqn.flux();
        }
    }

압력수정방정식을 풀고 있다. "fvc:div(phi)"는 아래와 같이 되어 있다.

$FOAM_SRC/finiteVolume/finiteVolume/fvc/fvcDiv.C

template<class Type>
tmp>GeometricField<Type, fvPatchField, volMesh> >
div
(
    const GeometricField& ssf
)
{
    return tmp >
    (
        new GeometricField
        (
            "div("+ssf.name()+')',
            fvc::surfaceIntegrate(ssf)
        )
    );
}

즉, "fvc::div(phi)"는 phi의 divergence가 아닌, phi를 이용해 divergence를 계산한다는 것을 의미한다. 위에서 기술된 "fvc::surfaceIntegrate()"를 쫒아가 보자.

finiteVolume/finiteVolume/fvc/fvcSurfaceIntegrate.C

template<class Type>
tmp<GeometricField<Type, fvPatchField, volMesh> >
surfaceIntegrate
(
    const GeometricField<Type, fvsPatfchField, surfaceMesh>& ssf
)
{
    const fvMesh& mesh = ssf.mesh();

    tmp <GeometricField<Type, fvPatchField, volMesh> > tvf
    (
        new GeometricField<Type, fvPatchField, volMesh>
        (
            IOobject
            (
                "surfaceIntegrate("+ssf.name()+')',
                ssf.instance(),
                mesh,
                IOobject::NO_READ,
                IOobject::NO_WRITE
            ),
            mesh,
            dimensioned
            (
                "0",
                ssf.dimensions()/dimVol,
                pTraits::zero
            ),
            zeroGradientFvPatchField<Type>::typeName
        )
    );
    GeometricField<Type, fvPatchField, volMesh>& vf = tvf();

    surfaceIntegrate(vf.internalField(), ssf);
    vf.correctBoundaryConditions();

    return tvf;
}

여기서 아래를 불러온다.

template<class Type>
void surfaceIntegrate
(
    Field<Type>& ivf,
    const GeometricField<Type, fvsPatchField, surfaceMesh>& ssf
)
{
    const fvMesh& mesh = ssf.mesh();

    const labelUList& owner = mesh.owner();
    const labelUList& neighbour = mesh.neighbour();

    const Field<Type>& issf = ssf;

    forAll(owner, facei)
    {
        ivf[owner[facei]] += issf[facei];
        ivf[neighbour[facei]] -= issf[facei];
    }

    forAll(mesh.boundary(), patchi)
    {
        const labelUList& pFaceCells =
            mesh.boundary()[patchi].faceCells();

        const fvsPatchField& pssf = ssf.boundaryField()[patchi];

        forAll(mesh.boundary()[patchi], facei)
        {
            ivf[pFaceCells[facei]] += pssf[facei];
        }
    }

    ivf /= mesh.V();
}

압력수정방정식으로 돌아가자.

        if (nonOrth == simple.nNonOrthCorr())
        {
            phi -= pEqn.flux();
        }

루프계산의 마지막에 유체 flux를 수정하고 있다.

압력수정방정식의 반복계산 루프는, 얼핏 같은 계산을 반복하는 의미없는 행위로 보일 지 모르지만, fvm::laplacian()이 음함수부분과 양함수부분을 둘 다 포함하므로, 매번 계산하는 식이 달라진다.

        p.relax();

p를 양함수로 relaxation하고 있다. SIMPLE법의 루프가 처음에 "p.storePrevIter()"를 불러오는 이유이다.

        U -= rAU*fvc::grad(p);
        U.correctBoundaryConditions();

유체 flux를 수정하고 경계조건에 적용시키고 있다.

simpleFoam.C로 돌아가자, 마지막부분에,

        turbulence->correct();

k-epsilon 모델의 경우를 보면,

$FOAM_SRC/turbulneceModels/incompressible/RAS?kEpsilon/kEpsilon.C

void kEpsilon::correct()
{
    RASModel::correct();

    if (!turbulence_)
    {
        return;
    }

    volScalarField G("RASModel::G", nut_*2*magSqr(symm(fvc::grad(U_))));

    // Update epsilon and G at the wall
    epsilon_.boundaryField().updateCoeffs();

    // Dissipation equation
    tmp epsEqn
    (
        fvm::ddt(epsilon_)
      + fvm::div(phi_, epsilon_)
      - fvm::Sp(fvc::div(phi_), epsilon_)
      - fvm::laplacian(DepsilonEff(), epsilon_)
     ==
        C1_*G*epsilon_/k_
      - fvm::Sp(C2_*epsilon_/k_, epsilon_)
    );

    epsEqn().relax();

    epsEqn().boundaryManipulate(epsilon_.boundaryField());

    solve(epsEqn);
    bound(epsilon_, epsilonMin_);


    // Turbulent kinetic energy equation
    tmp kEqn
    (
        fvm::ddt(k_)
      + fvm::div(phi_, k_)
      - fvm::Sp(fvc::div(phi_), k_)
      - fvm::laplacian(DkEff(), k_)
     ==
        G
      - fvm::Sp(epsilon_/k_, k_)
    );

    kEqn().relax();
    solve(kEqn);
    bound(k_, kMin_);


    // Re-calculate viscosity
    nut_ = Cmu_*sqr(k_)/epsilon_;
    nut_.correctBoundaryConditions();
}

ε 과 k 의 방정식을 풀어 난류점성을 계산한다.

위와 같이 살펴보아 "p.boundaryField(),updateCoeffs()"라던지 "U.correctBoundaryConditions()"의 역할은 잘 알 수 없었으나, kEpsilon::correct()를 보는한 방정식을 풀기 전에 "*.boundaryField().updateCoeffs()"를 불러와 U 의 보정이나 nut 를 계산하는 등의 직접변수를 계산할 때, 이후 "*.correctBoundaryCondition()"이 사용되는것을 확인할 수 있다. 

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