LoginSignup
4

More than 3 years have passed since last update.

posted at

updated at

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

使用環境

  • OpenFOAM v1806
  • Ubuntu 18.04

はじめに

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

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

解説

  • 変えた場所は"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();
        }
    }

    // 後略

最後に

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

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
What you can do with signing up
4