Help us understand the problem. What is going on with this article?

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

More than 1 year has passed since last update.

使用環境

  • 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();
        }
    }

    // 後略

最後に

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

Why do not you register as a user and use Qiita more conveniently?
  1. We will deliver articles that match you
    By following users and tags, you can catch up information on technical fields that you are interested in as a whole
  2. you can read useful information later efficiently
    By "stocking" the articles you like, you can search right away
Comments
Sign up for free and join this conversation.
If you already have a Qiita account
Why do not you register as a user and use Qiita more conveniently?
You need to log in to use this function. Qiita can be used more conveniently after logging in.
You seem to be reading articles frequently this month. Qiita can be used more conveniently after logging in.
  1. We will deliver articles that match you
    By following users and tags, you can catch up information on technical fields that you are interested in as a whole
  2. you can read useful information later efficiently
    By "stocking" the articles you like, you can search right away