chtMultiRegionFoamのカスタムソルバについて

[Satoshi]

お世話になっております。
中塚です。
現在、chtMultiRegionFoamで回転体の解析を行うため、カスタムソルばの作成を行っております。
simpleFoamとMRFsimpleFoamを見比べ、2つの違いを参考に下記のようにソル場を変更しました。
赤の太字の部分が追加した部分になります。
知識がないので、この程度の変更しかできないのですが、この様な変更でもソルバとして機能するのでしょうか？
他にも変更すべき場所などあればご指摘いただけると幸いです。
よろしくお願いします。


↓chtMultiRegionFoamカスタムソルバ


#include "fvCFD.H"
#include "basicRhoThermo.H"
#include "turbulenceModel.H"
#include "fixedGradientFvPatchFields.H"
#include "regionProperties.H"
#include "basicSolidThermo.H"
#include "radiationModel.H"
#include "MRFZones.H"


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

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

    regionProperties rp(runTime);

    #include "createFluidMeshes.H"
    #include "createSolidMeshes.H"

    #include "createFluidFields.H"
    #include "createSolidFields.H"

    #include "initContinuityErrs.H"
   
    MRFZones mrfZones(mesh);
    mrfZones.correctBoundaryVelocity(U);


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

        forAll(fluidRegions, i)
        {
            Info<< "\nSolving for fluid region "
                << fluidRegions[i].name() << endl;
            #include "setRegionFluidFields.H"
            #include "readFluidMultiRegionSIMPLEControls.H"
            #include "solveFluid.H"
        }

        forAll(solidRegions, i)
        {
            Info<< "\nSolving for solid region "
                << solidRegions[i].name() << endl;
            #include "setRegionSolidFields.H"
            #include "readSolidMultiRegionSIMPLEControls.H"
            #include "solveSolid.H"
        }

        runTime.write();

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

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

    return 0;
}


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

[ohbuchi]

pEqn.H, UEqn.Hも違う筈です。

On 6月29日, 午後11:47, Satoshi <nakatsukasatoshi05...@gmail.com> wrote:
> お世話になっております。
> 中塚です。
> 現在、chtMultiRegionFoamで回転体の解析を行うため、カスタムソルばの作成を行っております。
> simpleFoamとMRFsimpleFoamを見比べ、2つの違いを参考に下記のようにソル場を変更しました。
> 赤の太字の部分が追加した部分になります。
> 知識がないので、この程度の変更しかできないのですが、この様な変更でもソルバとして機能するのでしょうか？
> 他にも変更すべき場所などあればご指摘いただけると幸いです。
> よろしくお願いします。
>
> ↓chtMultiRegionFoamカスタムソルバ
>
> #include "fvCFD.H"
> #include "basicRhoThermo.H"
> #include "turbulenceModel.H"
> #include "fixedGradientFvPatchFields.H"
> #include "regionProperties.H"
> #include "basicSolidThermo.H"
> #include "radiationModel.H"

> *#include "MRFZones.H"*

>
> // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
> * //
>
> int main(int argc, char *argv[])
> {
> #include "setRootCase.H"
> #include "createTime.H"
>
> regionProperties rp(runTime);
>
> #include "createFluidMeshes.H"
> #include "createSolidMeshes.H"
>
> #include "createFluidFields.H"
> #include "createSolidFields.H"
>
> #include "initContinuityErrs.H"
>

> * MRFZones mrfZones(mesh);
> mrfZones.correctBoundaryVelocity(U);*

[Satoshi]

返信ありがとうございます。
UEqn.HについてはrhoPorousMRFSimpleFoam.Cを参考に
mrfZones.addCoriolis(UEqn());を加えてみました。

pEqn.Hについては mrfZones.relativeFlux(fvc::interpolate(***),*** ); のような一文を phi=
のあとに挿入していけばいいと思うのですが、よく分かりません。
アドバイスお願いします。


UEqn.H 
  // Solve the Momentum equation
    tmp<fvVectorMatrix> UEqn
    (
        fvm::div(phi, U)
      + turb.divDevRhoReff(U)
    );

    mrfZones.addCoriolis(UEqn());

    UEqn().relax();

    solve
    (
        UEqn()
     ==
        fvc::reconstruct
        (
            (
              - ghf*fvc::snGrad(rho)
              - fvc::snGrad(p_rgh)
            )*mesh.magSf()
        )
    );



pEqn.H

{
    rho = thermo.rho();
    rho = max(rho, rhoMin[i]);
    rho = min(rho, rhoMax[i]);
    rho.relax();

    volScalarField rAU(1.0/UEqn().A());
    surfaceScalarField rhorAUf("(rho*(1|A(U)))", fvc::interpolate(rho*rAU));

    p_rgh.boundaryField().updateCoeffs();
    U = rAU*UEqn().H();
    UEqn.clear();

    phi = fvc::interpolate(rho)*(fvc::interpolate(U) & mesh.Sf()); 
　mrfZones.relativeFlux(fvc::interpolate(rho),phi );　　　　　←ここに挿入?
    bool closedVolume = adjustPhi(phi, U, p_rgh);
    dimensionedScalar compressibility = fvc::domainIntegrate(psi);
    bool compressible = (compressibility.value() > SMALL);

    surfaceScalarField buoyancyPhi(rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf());
    phi -= buoyancyPhi;                                                      
　 mrfZones.relativeFlux(phi);　　　　　　　　　　 ←ここにも挿入？

    // Solve pressure
    for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
    {
        fvScalarMatrix p_rghEqn
        (
            fvm::laplacian(rhorAUf, p_rgh) == fvc::div(phi)
        );

        p_rghEqn.setReference
        (
            pRefCell,
            compressible ? getRefCellValue(p_rgh, pRefCell) : pRefValue
        );

        p_rghEqn.solve();

        if (nonOrth == nNonOrthCorr)
        {
            // Calculate the conservative fluxes
            phi -= p_rghEqn.flux()　　　　　　　　←この後にも？

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

            // Correct the momentum source with the pressure gradient flux
            // calculated from the relaxed pressure
            U -= rAU*fvc::reconstruct((buoyancyPhi + p_rghEqn.flux())/rhorAUf);
            U.correctBoundaryConditions();
        }
    }

    p = p_rgh + rho*gh;

    #include "continuityErrs.H"

    // For closed-volume cases adjust the pressure level
    // to obey overall mass continuity
    if (closedVolume && compressible)
    {
        p += (initialMass - fvc::domainIntegrate(thermo.rho()))
            /compressibility;
        p_rgh = p - rho*gh;
    }

    rho = thermo.rho();
    rho = max(rho, rhoMin[i]);
    rho = min(rho, rhoMax[i]);
    rho.relax();

    Info<< "Min/max rho:" << min(rho).value() << ' '
        << max(rho).value() << endl;

    // Update thermal conductivity
    kappa = thermo.Cp()*turb.alphaEff();
}



2012年6月30日土曜日 5時57分01秒 UTC+9 ohbuchi:

[Satoshi]

お世話になっております。

中塚です。

頂いた意見を参考に自分なりにソースファイルをいじり、

エラーががでずにコンパイルすることが出来ました。

これを使った計算などはまだしていませんが、

作成したソルバで間違ったところがあったら指摘していただけると助かります。

よろしくお願いいたします。

2012年6月30日土曜日 16時21分19秒 UTC+9 Satoshi:

[Satoshi]

お世話になっております。
中塚です。

chtMultiRegionFoamのカスタムソルバは完成し、それを用いて簡易的な解析を行っています。
作成したchtMultiRegionFoamのカスタムソルバについてはzipファイルでアップ致しました。
パスワードはfoamです。
http://www1.axfc.net/uploader/Ne/so/136324

計算が回ることには回るのですが、１００stepほどで止まってしまいます。
流体領域のfvSolutionのrelaxationFactorsもかなり小さくし、ある程度計算が回るようにはしているのですが、
根本的な解決にはなっておらずどうすればいいのか分からずにいます。
解決手段をご存知の方は語教授お願いいたします。

今回の解析結果の動画です。
http://www.youtube.com/watch?v=hZtCpl3GUj0&feature=youtu.be

流体領域のfvSolutionのrelaxationFactorsです。
relaxationFactors
{
     U 0.1;
     p 0.1;
     rho 0.1;
     h 0.1;
     k 0.1;
     omega 0.1;


2012年7月2日月曜日 21時44分11秒 UTC+9 Satoshi:

[Satoshi]

お世話になっております。
中塚です。

計算が止まる直前のタイムステップの表示は下記のようになっていました。


Time = 0.0084


Solving for fluid region case
DILUPBiCG:  Solving for Ux, Initial residual = 0.0042200332, Final residual = 9.4022068e-06, No Iterations 1
DILUPBiCG:  Solving for Uy, Initial residual = 0.0036325088, Final residual = 1.0635911e-05, No Iterations 1
DILUPBiCG:  Solving for Uz, Initial residual = 0.013924242, Final residual = 3.455069e-05, No Iterations 1
DILUPBiCG:  Solving for h, Initial residual = 0.02799126, Final residual = 4.3352043e-05, No Iterations 1
Min/max T:-15.559199 1378.0684
GAMG:  Solving for p_rgh, Initial residual = 0.033198146, Final residual = 0.0002088189, No Iterations 4
time step continuity errors : sum local = 1.1026978e-05, global = 5.307957e-07, cumulative = 6.3598185e-05
Min/max rho:0.99985666 1.0298567

Solving for solid region Solid
DICPCG:  Solving for T, Initial residual = 9.5590812e-07, Final residual = 9.5590812e-07, No Iterations 0
DICPCG:  Solving for T, Initial residual = 9.5614709e-07, Final residual = 9.5614709e-07, No Iterations 0
Min/max T:300 500
ExecutionTime = 73.2 s  ClockTime = 74 s

Time = 0.0085


Solving for fluid region case
DILUPBiCG:  Solving for Ux, Initial residual = 0.0043183715, Final residual = 9.4957365e-06, No Iterations 1
DILUPBiCG:  Solving for Uy, Initial residual = 0.0037662967, Final residual = 1.0961174e-05, No Iterations 1
DILUPBiCG:  Solving for Uz, Initial residual = 0.014021344, Final residual = 3.4541208e-05, No Iterations 1
DILUPBiCG:  Solving for h, Initial residual = 0.028325685, Final residual = 4.2171554e-05, No Iterations 1


--> FOAM FATAL ERROR:
Maximum number of iterations exceeded

    From function specieThermo<Thermo>::T(scalar f, scalar T0, scalar (specieThermo<Thermo>::*F)(const scalar) const, scalar (specieThermo<Thermo>::*dFdT)(const scalar) const) const
    in file /home/opencfd/OpenFOAM/OpenFOAM-2.1.1/src/thermophysicalModels/specie/lnInclude/specieThermoI.H at line 69.

FOAM aborting

#0  Foam::error::printStack(Foam::Ostream&) in "/opt/openfoam211/platforms/linux64GccDPOpt/lib/libOpenFOAM.so"
#1  Foam::error::abort() in "/opt/openfoam211/platforms/linux64GccDPOpt/lib/libOpenFOAM.so"
#2  Foam::specieThermo<Foam::hConstThermo<Foam::perfectGas> >::T(double, double, double (Foam::specieThermo<Foam::hConstThermo<Foam::perfectGas> >::*)(double) const, double (Foam::specieThermo<Foam::hConstThermo<Foam::perfectGas> >::*)(double) const, double (Foam::specieThermo<Foam::hConstThermo<Foam::perfectGas> >::*)(double) const) const in "/opt/openfoam211/platforms/linux64GccDPOpt/lib/libbasicThermophysicalModels.so"
#3  Foam::hRhoThermo<Foam::pureMixture<Foam::constTransport<Foam::specieThermo<Foam::hConstThermo<Foam::perfectGas> > > > >::calculate() in "/opt/openfoam211/platforms/linux64GccDPOpt/lib/libbasicThermophysicalModels.so"
#4  Foam::hRhoThermo<Foam::pureMixture<Foam::constTransport<Foam::specieThermo<Foam::hConstThermo<Foam::perfectGas> > > > >::correct() in "/opt/openfoam211/platforms/linux64GccDPOpt/lib/libbasicThermophysicalModels.so"
#5 
 in "/home/taro/OpenFOAM/taro-2.1.1/platforms/linux64GccDPOpt/bin/chtMRFFoam"
#6  __libc_start_main in "/lib/x86_64-linux-gnu/libc.so.6"
#7 
 in "/home/taro/OpenFOAM/taro-2.1.1/platforms/linux64GccDPOpt/bin/chtMRFFoam"
中止 (コアダンプ)


2012年7月18日水曜日 18時49分09秒 UTC+9 Satoshi: