OpenFOAMで場所ごとの残差を出力

使用環境

• OpenFOAM v1806
• Ubuntu 18.04

はじめに

　OpenFOAMで定常計算をしていてなかなか収束しないときに、どの場所で発散しているかを知りたいことがあるかもしれない。そんなときはfvMatrix::residual()という関数を使うとその時の各セルの残差を返してくれる。
　以下はsimpleFoamに追加する場合の例。

例

simpleFoamのpとUの初期残差(initial residual)を各timeStepで出力する。
https://github.com/inabower/residualSimpleFoam

解説

• 変えた場所は"createFields.H"と"(p|U)Eqn.H"。
• fvMatrix::residual()という関数でその時の各セルにおける左右辺の差（残差）を求める。
• vol(Scalar|Vector)Fieldに無次元のfieldを代入する際にはprimitiveFieldRef()関数を使用する。
• 今回のコードはinitial residualの出力ですが、各Eqn.Hの中のsolve()の後に置けばfinal residualを求めることも可能。

residualSimpleFoam/createFields.H

// 元々のcreateFields.Hに以下を追記

Info<< "Creating residual field of p and U\n" << endl;
volScalarField rp // pの残差
(
    IOobject
    (
        "residual_p",
        runTime.timeName(),
        mesh,
        IOobject::NO_READ,
        IOobject::AUTO_WRITE
    ),
    mesh,
    dimensionedScalar("resp", dimless, 0.)
);

volVectorField rU // Uの残差
(
    IOobject
    (
        "residual_U",
        runTime.timeName(),
        mesh,
        IOobject::NO_READ,
        IOobject::AUTO_WRITE
    ),
    mesh,
    dimensionedVector("resU", dimless, Zero)
);

residualSimpleFoam/UEqn.H

    // Momentum predictor

    MRF.correctBoundaryVelocity(U);

    tmp<fvVectorMatrix> tUEqn
    (
        fvm::div(phi, U)
      + MRF.DDt(U)
      + turbulence->divDevReff(U)
     ==
        fvOptions(U)
    );
    fvVectorMatrix& UEqn = tUEqn.ref();
    
    // ここを追加（最終残差を見たい場合はsolveの後に置く）
    rU.primitiveFieldRef() = UEqn.residual();

    UEqn.relax();

    fvOptions.constrain(UEqn);

    if (simple.momentumPredictor())
    {
        solve(UEqn == -fvc::grad(p));

        fvOptions.correct(U);
    }

residualSimpleFoam/pEqn.H

    // 前略

    tUEqn.clear();

    // ここを追加（初期残差が欲しい場合）
    bool isInitialResidual = true;
    
    // Non-orthogonal pressure corrector loop
    while (simple.correctNonOrthogonal())
    {
        fvScalarMatrix pEqn
        (
            fvm::laplacian(rAtU(), p) == fvc::div(phiHbyA)
        );

        // ここを追加（初期残差が欲しい場合）
        if(isInitialResidual)
        {
            rp.primitiveFieldRef() = pEqn.residual();
            isInitialResidual = false;
        }

        pEqn.setReference(pRefCell, pRefValue);

        pEqn.solve();

        if (simple.finalNonOrthogonalIter())
        {
            phi = phiHbyA - pEqn.flux();
        }
    }

    // 後略

最後に

　全然収束しないときに役立つかもしれません。