水の自然対流の解析について

[hiro]

こんにちは，OpenFOAM初心者の学生hiroです．

OpenFOAM2.3.0を使い始めて，約1カ月です．

今，行き詰っておりますので，どなたかご教授ください．

自分のやりたい解析としては，水の自然対流を解析したいと思っております．

しかし，OpenFOAM2.3.0では似たようなチュートリアルが無かったため，

古いバーション（1.6）のチュートリアルからコピペしたbuoyantSimpleFoamの「hotRoom」を

使って自然対流の解析を行いました．

しかし，エラーが出てしまい，行き詰っております．

このエラーの原因が　①バージョンが違うから

　　　　　　　　　　　　　　②コピー＆ペーストしたから

　　　　　　　　　　　　　　③その他　　　　　　　　　　　　…なのかがわかりません．

エラーの原因がわかる方がいれば，アドバイスお願いします．

一応，エラーの文章も記載しておきますので，何卒よろしくお願いします．

--> FOAM FATAL ERROR:
Unknown rhoThermo type hPsiThermo<pureMixture<constTransport<specieThermo<hConstThermo<perfectGas>>>>>

Valid rhoThermo types are:

52
(
heRhoThermo<homogeneousMixture<const<hConst<incompressiblePerfectGas<specie>>,sensibleEnthalpy>>>
heRhoThermo<homogeneousMixture<const<hConst<perfectGas<specie>>,sensibleEnthalpy>>>
heRhoThermo<homogeneousMixture<sutherland<janaf<incompressiblePerfectGas<specie>>,sensibleEnthalpy>>>
heRhoThermo<homogeneousMixture<sutherland<janaf<perfectGas<specie>>,sensibleEnthalpy>>>
heRhoThermo<inhomogeneousMixture<const<hConst<incompressiblePerfectGas<specie>>,sensibleEnthalpy>>>
heRhoThermo<inhomogeneousMixture<const<hConst<perfectGas<specie>>,sensibleEnthalpy>>>
heRhoThermo<inhomogeneousMixture<sutherland<janaf<incompressiblePerfectGas<specie>>,sensibleEnthalpy>>>
heRhoThermo<inhomogeneousMixture<sutherland<janaf<perfectGas<specie>>,sensibleEnthalpy>>>
heRhoThermo<multiComponentMixture<const<hConst<incompressiblePerfectGas<specie>>,sensibleEnthalpy>>>
heRhoThermo<multiComponentMixture<const<hConst<incompressiblePerfectGas<specie>>,sensibleInternalEnergy>>>
heRhoThermo<multiComponentMixture<const<hConst<perfectGas<specie>>,sensibleEnthalpy>>>
heRhoThermo<multiComponentMixture<const<hConst<perfectGas<specie>>,sensibleInternalEnergy>>>
heRhoThermo<multiComponentMixture<polynomial<hPolynomial<icoPolynomial<specie>>,sensibleEnthalpy>>>
heRhoThermo<multiComponentMixture<polynomial<hPolynomial<icoPolynomial<specie>>,sensibleInternalEnergy>>>
heRhoThermo<multiComponentMixture<sutherland<janaf<incompressiblePerfectGas<specie>>,sensibleEnthalpy>>>

…)

     From function rhoThermo::New

     in file lnInclude/basicThermoTemplates.C at line 123.

FOAM exiting

以上になります．よろしくお願いします．

[ohbuchi]

こんにちは。
OpenFOAM-2.2から熱物性モデルの記述方法が変わりました。
それが原因で発生したエラーです。
buoyantPimpleFoam/hotRoomチュートリアルをbuoyantSimpleFoamで計算すれば
お望みの計算が簡単に実行できると思います。


2014年7月23日水曜日 19時16分05秒 UTC+9 hiro:

[hiro]

こんにちは．

ohbuchiさん，お返事ありがとうございます．

OpenFOAM初心者のhiroです．

ohbuchiさんのアドバイス通り

buoyantPimpleFoam/hotRoomのチュートリアルを

buoyantSimpleFoamのフォルダにコピペして計算させました．

しかし，以下のエラーが出てきてしまいました．

--> FOAM FATAL IO ERROR:
keyword SIMPLE is undefined in dictionary "/root/OpenFOAM/root-2.3.0/run/tutorials/heatTransfer/buoyantSimpleFoam/hotRoom2/system/fvSolution"

file: /root/OpenFOAM/root-2.3.0/run/tutorials/heatTransfer/buoyantSimpleFoam/hotRoom2/system/fvSolution from line 22 to line 62.

    From function dictionary::subDict(const word& keyword) const
    in file db/dictionary/dictionary.C at line 643.

FOAM exiting

以上になります．

このエラーの原因・意味はどういうことなのでしょうか？

またbuoyantPimpleFoam/hotRoomチュートリアルをbuoyantSimpleFoamで

計算させて，定常（Simple）状態の結果が得られるのでしょうか？

それとも私がohbuchiさんのアドバイスを間違って解釈してしまったのでしょうか？

初心者なので，どうかご教授ください．

もしトゲのある言い方になっていたら，申し訳ございません．

[ohbuchi]

エラーメッセージを良く読んでください。
ｆｖSolutionファイルにキーワードSIMPLEの定義がないと書いてあります。
buoyantSimpleFoamのチュートリアル中のｆｖSolutionファイルから、下記の記述を
コピーすればOKです。
SIMPLE
{
    nNonOrthogonalCorrectors 0;
    pRefCell        0;
    pRefValue       0;

    residualControl
    {
        p_rgh           1e-2;
        U               1e-3;
        h               1e-3;

        // possibly check turbulence fields
        "(k|epsilon|omega)" 1e-3;
    }
}

relaxationFactors
{
    fields
    {
        rho             1.0;
        p_rgh           0.7;
    }
    equations
    {
        U               0.2;
        h               0.2;
        "(k|epsilon|R)" 0.5;
    }
}


2014年7月24日木曜日 12時28分44秒 UTC+9 hiro:

[hiro]

ohbuchiさん，ありがとうございました．

言われた通り，実行したら計算ができました．

これからはしっかり，エラーを読んでいきたいと思います．

今後も何かと質問をすると思いますが，よろしくお願いします．

[hiro]

流体を水に設定にするには，どうしたらよいのでしょうか？

thermophysicalPropertiesを変更するのはわかるんですけど，

物性値以外，どのように変更したいいのか分かりません．

他のtutorialを確認しても，ほとんど空気ばかりで

空気と水を用いているtutorialでも水は”気体”として扱われているので

”液体”として扱うには，どうしたらよろしいでしょうか？

[ohbuchi]

chtMultiRegionFoamのチュートリアルでbotomWaterの部分のthermophysicalPropertiesを参考にすると
密度一定の場合の液体解析ができます。しかし、密度一定で浮力を扱うにはbuoyantBoussinesqSimpleFoamに
する必要があると思います。
液体で、圧縮性を考慮するには、compressibleInterFoamのチュートリアルであるdepthCharge2DチュートリアルのthermophysicalProperties.waterが参考になります。


2014年8月5日火曜日 14時35分24秒 UTC+9 hiro:

[hiro]

ohbuchi様

アドバイスありがとうございます．

流体の勉強を含め， buoyantBoussinesqSimpleFoamでやり直しました．

水の物性値をtransportPropertiesに書き込み，計算を実行できました．

ありがとうございます．

しかし，計算途中でエラーになっていまいました．

そこでblockMeshDictをcheckMeshで確認したところ（後半部にて）

Checking geometry...

    Overall domain bounding box (0 0 0) (0.08 0.08 0.08)

    Mesh (non-empty, non-wedge) directions (1 1 1)

    Mesh (non-empty) directions (1 1 1)

    Boundary openness (-3.41079e-17 -7.07985e-17 5.32252e-17) OK.

    Max cell openness = 2.04058e-16 OK.

    Max aspect ratio = 9.29557 OK.

    Minimum face area = 1.91144e-06. Maximum face area = 1.61179e-05.  Face area magnitudes OK.

    Min volume = 7.64576e-09. Max volume = 6.44717e-08.  Total volume = 0.000402298.  Cell volumes OK.

    Mesh non-orthogonality Max: 73.4401 average: 18.7088

   *Number of severely non-orthogonal (> 70 degrees) faces: 160.

    Non-orthogonality check OK.

  <<Writing 160 non-orthogonal faces to set nonOrthoFaces

    Face pyramids OK.

    Max skewness = 0.994343 OK.

    Coupled point location match (average 0) OK.

 

Mesh OK.

…とメッシュの非直行性に問題（メッシュの歪みが大きい）と出たことで

ソルバーの計算を継続できなかったとわかりました．

しかし，今の作成モデルをどのように改良したら，低い非直行性にできるか分かりません．

作成したblockMeshDictを記載しますので，助言を頂けませんか？

ちなみに自力で作った円筒モデルなので分かりにくいかもしれません…

convertToMeters 0.001;

vertices
(
    (40 0 0) //0
    (68.3 0 11.7) //1
    (80 0 40) //2
    (68.3 0 68.3) //3
    (40 0 80) //4
    (11.7 0 68.3) //5
    (0 0 40) //6
    (11.7 0 11.7) //7
    (40 0 40) //8
    (40 80 0) //9
    (68.3 80 11.7) //10
    (80 80 40) //11
    (68.3 80 68.3) //12
    (40 80 80) //13
    (11.7 80 68.3) //14
    (0 80 40) //15
    (11.7 80 11.7) //16
    (40 80 40) //17
);

edges
(
    arc 7 0 (24.7 0 3)
    arc 0 1 (55.3 0 3)
    arc 1 2 (77 0 24.7)
    arc 2 3 (77 0 55.3)
    arc 3 4 (55.3 0 77)
    arc 4 5 (24.7 0 77)
    arc 5 6 (3 0 55.3)
    arc 6 7 (3 0 24.7)
    arc 16 9 (24.7 80 3)
    arc 9 10 (55.3 80 3)
    arc 10 11 (77 80 24.7)
    arc 11 12 (77 80 55.3)
    arc 12 13 (55.3 80 77)
    arc 13 14 (24.7 80 77)
    arc 14 15 (3 80 55.3)
    arc 15 16 (3 80 24.7)
);

blocks
(
    hex (7 0 9 16 6 8 17 15) (10 20 10) simpleGrading (1 1 1)
    hex (0 1 10 9 8 2 11 17) (10 20 10) simpleGrading (1 1 1)
    hex (6 8 17 15 5 4 13 14) (10 20 10) simpleGrading (1 1 1)
    hex (8 2 11 17 4 3 12 13) (10 20 10) simpleGrading (1 1 1)
);

boundary
(
    fixedWalls
    {
        type wall;
        faces
        (
            (0 9 10 1)
            (0 9 16 7)
            (1 10 11 2)
            (7 16 15 6)
            (2 11 12 3)
            (6 15 14 5)
            (3 12 13 4)
            (5 14 13 4)
        );
    }

    hot
    {
        type wall;
        faces
        (
            (0 1 2 8)
            (0 7 6 8)
            (2 8 4 3)
            (8 6 5 4)
        );
    }

    cold
    {
        type wall;
        faces
        (
            (9 10 11 17)
            (9 16 15 17)
            (11 17 13 12)
            (17 15 14 13)
        );
    }
);

mergePatchPairs
(
);

アドバイスのほう，よろしくお願いします．

[ohbuchi]

こんにちは。
円柱を１つのブロックで表現しようとすると、ブロックのコーナー角が大きくなり、どうしても
非直交性が悪化します。
円柱の中心部に角柱を設け、周囲に扇形状のブロックを４つ配する様な構造にすると直交性が
大幅に改善されます。

E.Moguraさんが公開しているマクロを使うと簡単に作成することができます。
http://mogura7.zenno.info/~et/xoops/modules/wordpress/index.php?p=408

ご参考まで。


2014年8月8日金曜日 3時02分26秒 UTC+9 hiro:

[hiro]

こんにちは．

先日，ohbuchi様から助言頂いたとおり

E.Moguraさんが公開しているマクロを使った円筒モデルを用いることで

blockMeshやcheckMeshが問題なく，通るようになりました．

しかし，また新たな問題点が出てきました．

計算途中でエラーが表示され，それ以上の計算がされなくなってしまいました．

原因を調べ，fvSchemsかfvSolutionの中を変更すればできるという考えていますが，

どこをどのように変更すればいいか分かりません…

そなたか，ご教授お願いします．

下記にエラー内容を記載しますので，よろしくお願いします．

Time = 2530

 

DILUPBiCG:  Solving for Ux, Initial residual = 0.0772779, Final residual = 7.10049e-05, No Iterations 1

DILUPBiCG:  Solving for Uy, Initial residual = 0.0973584, Final residual = 0.000137618, No Iterations 1

DILUPBiCG:  Solving for Uz, Initial residual = 0.101135, Final residual = 0.000623943, No Iterations 1

DILUPBiCG:  Solving for T, Initial residual = 0.0153583, Final residual = 0.000420818, No Iterations 1

DICPCG:  Solving for p_rgh, Initial residual = 0.557855, Final residual = 0.00459961, No Iterations 76

time step continuity errors : sum local = 4.38405e+90, global = -4.33121e+90, cumulative = -1.41442e+91

#0  Foam::error::printStack(Foam::Ostream&) in "/opt/openfoam230/platforms/linuxGccDPOpt/lib/libOpenFOAM.so"

#1  Foam::sigFpe::sigHandler(int) in "/opt/openfoam230/platforms/linuxGccDPOpt/lib/libOpenFOAM.so"

#2  Uninterpreted:

#3  Foam::PBiCG::solve(Foam::Field<double>&, Foam::Field<double> const&, unsigned char) const in "/opt/openfoam230/platforms/linuxGccDPOpt/lib/libOpenFOAM.so"

#4  Foam::fvMatrix<double>::solveSegregated(Foam::dictionary const&) in "/opt/openfoam230/platforms/linuxGccDPOpt/lib/libfiniteVolume.so"

#5  Foam::fvMatrix<double>::solve(Foam::dictionary const&) in "/opt/openfoam230/platforms/linuxGccDPOpt/bin/buoyantBoussinesqSimpleFoam"

#6  Foam::fvMatrix<double>::solve() in "/opt/openfoam230/platforms/linuxGccDPOpt/lib/libincompressibleRASModels.so"

#7  Foam::incompressible::RASModels::kEpsilon::correct() in "/opt/openfoam230/platforms/linuxGccDPOpt/lib/libincompressibleRASModels.so"

#8 

 in "/opt/openfoam230/platforms/linuxGccDPOpt/bin/buoyantBoussinesqSimpleFoam"

#9  __libc_start_main in "/lib/i386-linux-gnu/libc.so.6"

#10 

 in "/opt/openfoam230/platforms/linuxGccDPOpt/bin/buoyantBoussinesqSimpleFoam"

浮動小数点例外 (コアダンプ)

2014年8月8日金曜日 8時53分03秒 UTC+9 ohbuchi:

[Youhei Takagi]

hiro様

高木と申します。

continuity errorが大きくなっているので発散しているのだと思いますが、
これだけの情報では問題点がわかりません。格子解像度、時間刻み、物性値、
初期条件、境界条件等すべて提供されないと妥当な計算を行っているか
どうか判断できません。

ご存知だと思いますが、空気のプラントル数は0.7、水の場合は7ぐらいですので、
流体を変更したら計算格子の品質を変えるのが普通だと思います。

以上ご参考までに。

2014年9月4日 9:01 hiro <mg.ra...@gmail.com>:

--
このメールは Google グループのグループ「OpenFOAM」に登録しているユーザーに送られています。
このグループから退会し、グループからのメールの配信を停止するには openfoam+u...@googlegroups.com にメールを送信してください。
このグループに投稿するには open...@googlegroups.com にメールを送信してください。
http://groups.google.com/group/openfoam からこのグループにアクセスしてください。
その他のオプションについては https://groups.google.com/d/optout にアクセスしてください。

[hiro]

高木様，アドバイスありがとうございます.

ちなみに計算格子の品質を変えるのが普通とのことですが，
この計算格子とはblockMeshを変更すれば，よろしいのでしょうか？

m4を用いて作成した円柱モデルなので，直交性は大幅に改善されたと
思ったのですが…


2014年9月4日木曜日 9時26分22秒 UTC+9 yotakagi:

[Youhei Takagi]

hiro様

ご存知だと思いますが、プラントル数は速度境界層と温度境界層の
厚みの比率を表しますので、高プラントル数流体は低プラントル数流体に
比べて細かい計算格子が必要となります。

ただし、今回の発散の原因は計算格子が原因とは限らないので
(前述のとおりすべての情報が提示されていないので)、まずはご自身で
かならず計算が走る条件で検証されるのがよいかと思います。

2014年9月4日 18:18 hiro <mg.ra...@gmail.com>:

[hiro]

yotakagi 様，ありがとうございます．

まずは計算格子の大きさを変えてやってみます．

調べてみると，fvSchemsで指定されているスキームを

別のものに変える方法もあるみたいですね．

他にも，境界条件等も確認して，試したいと思います．

改めて，アドバイスありがとうございました．

2014年9月4日木曜日 19時48分01秒 UTC+9 yotakagi:

[hiro]

皆様

こんにちは，いつもお世話になっております．

前回の計算エラーの原因が，初期条件の設定に問題があることが分かり

設定を変えたら，計算エラーが出ずに計算を終えることができました．

ohbuchi様，高木様 ありがとうございました．

しかし，計算結果が定常になるまで

かなり時間・容量がかかってしまうことが分かりました．

そこで今，円柱で行っている計算を

半円柱に変え，symmetryPlane (対称面)を使い

円柱と仮定して計算を行おうと考えています．

そこで一つ質問があります．

モデルが対称であるのと同様に

自然対流の計算結果も対称的になってしまうのでしょうか？

つまり，symmetryPlaneを使った半円柱モデルの計算結果では

自然対流の仕方が鏡で見るような対称的になってしまうのかと

いう質問です．

わかりづらい質問で誠に申し訳ありませんが，

よろしくお願いします．

2014年9月5日金曜日 13時52分27秒 UTC+9 hiro:

[ohbuchi]

こんにちは。
SymmetryPlane境界では境界に垂直な方向の勾配を０とするため、非対称な解は得られません。
非対称な結果が想定される場合には、妥当性のない境界条件設定となります。


2014年9月9日火曜日 8時18分26秒 UTC+9 hiro:

[hiro]

こんにちは．

ohbuchi様，いつも貴重なアドバイスありがとうございます．

大変勉強になりました．もう少し考えていきたいと思います．

一応，三次元の半円柱モデルで自然対流の数値シミュレーションを行ったのですが，

Time=20000まで計算を行っても，定常解を得られることができませんでした．

そこで簡単に

二次元のモデルで計算を行いました．

しかし，二次元にしても定常解を得ることができませんでした．

（同様にTime=20000で行いました）

二次元にしても，定常解が得られないのは

納得いきませんが，どの部分が問題なのかわかりません．

そこで初期条件などを記載していますので，

ご指摘，お願いできないでしょうか？

------blockMeshDict-------

convertToMeters 0.001;

vertices
(
    (0 0 0)
    (80 0 0)
    (80 1 0)
    (0 1 0)
    (0 0 50)
    (80 0 50)
    (80 1 50)
    (0 1 50)

);

edges
(
);

blocks
(
    hex (0 1 2 3 4 5 6 7) (40 1 20) simpleGrading (1 1 1)
);

boundary
(
    fixedWalls
    {
        type wall;
        faces
        (

            (0 1 5 4)
            (1 2 6 5)
            (0 3 7 4)
            (3 2 6 7)
        );
    }

    hot
    {
        type wall;
        faces
        (

            (0 1 2 3)
        );
    }

    cold
    {
        type wall;
        faces
        (

            (4 5 6 7)

        );
    }
);

mergePatchPairs
(
);

------controlDict------

application     buoyantBoussinesqSimpleFoam;

startFrom       latestTime;

startTime       0;

stopAt          endTime;

endTime         20000;

deltaT          1;

writeControl    timeStep;

writeInterval   50;

purgeWrite      0;

writeFormat     ascii;

writePrecision  6;

writeCompression off;

timeFormat      general;

timePrecision   6;

runTimeModifiable true;

------fvSchemes------

ddtSchemes
{
    default         steadyState;
}

gradSchemes
{
    default         Gauss linear;
}

divSchemes
{
    default         none;

    div(phi,U)      bounded Gauss upwind;
    div(phi,T)      bounded Gauss upwind;
    div(phi,k)      bounded Gauss upwind;
    div(phi,epsilon) bounded Gauss upwind;
    div((nuEff*dev(T(grad(U))))) Gauss linear;
}

laplacianSchemes
{
    default         Gauss linear corrected;
}

interpolationSchemes
{
    default         linear;
}

snGradSchemes
{
    default         corrected;
}

fluxRequired
{
    default         no;
    p_rgh           ;
}

------fvSolution------
solvers
{
    p_rgh
    {
        solver          PCG;
        preconditioner  DIC;
        tolerance       1e-08;
        relTol          0.01;
    }

    "(U|T|k|epsilon|R)"
    {
        solver          PBiCG;
        preconditioner  DILU;
        tolerance       1e-05;
        relTol          0.1;
    }
}

SIMPLE
{
    nNonOrthogonalCorrectors 0;
    pRefCell        0;
    pRefValue       0;

    residualControl
    {
        p_rgh           1e-2;
        U               1e-4;
        T               1e-2;

        // possibly check turbulence fields
        "(k|epsilon|omega)" 1e-3;
    }
}

relaxationFactors
{
    fields
    {
        p_rgh           0.7;
    }
    equations
    {
        U               0.3;
        T               1;
        "(k|epsilon|R)" 0.7;
    }
}

------setFieldsDict------

defaultFieldValues
(
    volScalarFieldValue T 300
);

regions
(
    // Set patch values (using ==)
    boxToFace
    {
        box (4.5 -1000 4.5) (5.5 1e-5 5.5);

        fieldValues
        (
            volScalarFieldValue T 600
        );
    }
);

------g------

dimensions      [0 1 -2 0 0 0 0];
value           ( 0 0 -9.81 );

------RASProperties------

RASModel        laminar;

turbulence      on;

printCoeffs     on;

------thermophysicalProperties------

thermoType
{
    type            heRhoThermo;
    mixture         pureMixture;
    transport       const;
    thermo          hConst;
    equationOfState rhoConst;
    specie          specie;
    energy          sensibleEnthalpy;
}

mixture
{
    specie
    {
        nMoles          1;
        molWeight       18;
    }
    equationOfState
    {
        rho             1000;
    }
    thermodynamics
    {
        Cp              4181;
        Hf              0;
    }
    transport
    {
        mu              959e-6;
        Pr              6.62;
    }
}

------transportProperties------

transportModel Newtonian;

// Laminar viscosity
nu              nu [0 2 -1 0 0 0 0] 8.71e-07;

// Thermal expansion coefficient
beta            beta [0 0 0 -1 0 0 0] 2.57e-04;

// Reference temperature
TRef            TRef [0 0 0 1 0 0 0] 293;

// Laminar Prandtl number
Pr              Pr [0 0 0 0 0 0 0] 5.77;

// Turbulent Prandtl number
Prt             Prt [0 0 0 0 0 0 0] 1;

------0/T------

dimensions      [0 0 0 1 0 0 0];

internalField   uniform 293;

boundaryField
{
    fixedWalls
    {
        type            zeroGradient;
    }

    hot
    {
        type            fixedValue;
        value           uniform 313.15; // 40 degC
    }
    cold
    {
        type            fixedValue;
        value           uniform 288.15; // 15 degC
    }
}

------0/U------

dimensions      [0 1 -1 0 0 0 0];

internalField   uniform (0 0 0);

boundaryField
{
    fixedWalls
    {
        type            fixedValue;
        value           uniform (0 0 0);
    }

    hot
    {
        type            fixedValue;
        value           uniform (0 0 0);
    }

    cold
    {
        type            fixedValue;
        value           uniform (0 0 0);
    }
}

------0/alphat------

dimensions      [0 2 -1 0 0 0 0];

internalField   uniform 0;

boundaryField
{
    fixedWalls
    {
        type            alphatJayatillekeWallFunction;
        Prt             0.85;
        value            $internalField;
    }
    hot
    {
        type            alphatJayatillekeWallFunction;
        Prt             0.85;
        value           $internalField;
    }
    cold
    {
        type            alphatJayatillekeWallFunction;
        Prt             0.85;
        value           $internalField;
    }
}

------0/epsilon------

dimensions      [0 2 -3 0 0 0 0];

internalField   uniform 0.01;

boundaryField
{
    fixedWalls
    {
        type            epsilonWallFunction;
        value            $internalField;
    }
    hot
    {
        type            epsilonWallFunction;
        value            $internalField;
    }
    cold
    {
        type            epsilonWallFunction;
        value            $internalField;
    }
}

------0/k------

dimensions      [0 2 -2 0 0 0 0];

internalField   uniform 0.1;

boundaryField
{
    fixedWalls
    {
        type            kqRWallFunction;
        value            $internalField;
    }
    hot
    {
        type            kqRWallFunction;
        value            $internalField;
    }
    cold
    {
        type            kqRWallFunction;
        value            $internalField;
    }
}

------0/nut------

dimensions      [0 2 -1 0 0 0 0];

internalField   uniform 0;

boundaryField
{
    fixedWalls
    {
        type            nutkWallFunction;
        value           uniform 0;
    }
    hot
    {
        type            nutkWallFunction;
        value           uniform 0;
    }
    cold
    {
        type            nutkWallFunction;
        value           uniform 0;
    }
}

------0/p------

dimensions      [0 2 -2 0 0 0 0];

internalField   uniform 0;

boundaryField
{
    fixedWalls
    {
        type            calculated;
        value           $internalField;
    }

    hot
    {
        type            calculated;
        value           $internalField;
    }

    cold
    {
        type            calculated;
        value           $internalField;
    }
}

------0/p_rgh------

dimensions      [0 2 -2 0 0 0 0];

internalField   uniform 0;

boundaryField
{
    fixedWalls
    {
        type            fixedFluxPressure;
        rho             rhok;
        value           $internalField;
    }

    hot
    {
        type            fixedFluxPressure;
        rho             rhok;
        value           $internalField;
    }

    cold
    {
        type            fixedFluxPressure;
        rho             rhok;
        value           $internalField;
    }
}

以上が

自分の作成した水の自然対流シミュレーション（二次元）の初期条件等になります．

どうぞ，よろしくおねがいします．

2014年9月9日火曜日 8時58分06秒 UTC+9 ohbuchi:

[ohbuchi]

まず、これでは2次元計算になっていません。
２次元計算では鉛直軸をｙ軸にとり、z軸に垂直な面はempty型のパッチとします。
blockMeshDictは下記の様になります。

convertToMeters 0.001;
vertices
(
    (0  0  0)
    (80 0  0)

    (80 50 0)
    (0  50 0)
    (0  0  1)
    (80 0  1)
    (80 50 1)
    (0  50 1)

);
edges
(
);
blocks
(

    hex (0 1 2 3 4 5 6 7) (40 20 1) simpleGrading (1 1 1)

);
boundary
(
    fixedWalls
    {
        type wall;
        faces
        (

            (0 4 7 3)
            (1 2 6 5)

        );
    }
    hot
    {
        type wall;
        faces
        (

            (0 1 5 4)

        );
    }
    cold
    {
        type wall;
        faces
        (

            (3 7 6 2)
        );
    }
    empty
    {
        type  empty;
        faces
        (
            (0 3 2 1)

            (4 5 6 7)
        );
    }
);
mergePatchPairs
(
);

関連してgは（0 -9.81 0)にします。
また、０ディレクトリの全ファイルにemptyパッチを追加します。



2014年9月11日木曜日 12時14分51秒 UTC+9 hiro:

[hiro]

ohbuchi様，ありがとうございます．

アドバイス通り，blockMeshを変更したら定常解を得られました．（キャプチャ8）

しかし，実際の現象と比較するとかなり違うものになっています．

むしろ定常解を得られなかったときのほうが

実際の現象に近い結果が得られていました．（キャプチャ9）

今回得られた定常解の結果を見てみると

右側側面に温度の高い流体が沿って流れていると分かります．

この結果からだと側面が断熱されていない可能性が考えれます．

初期設定では，上面は冷却（15℃），下面は加熱（40℃），そして側面は断熱されているはずです．

そこで質問なのですが，

初期設定の温度設定（0/T）は，始めのときだけで

それ以外ではリセットされているのでしょうか？

つまり，0/Tで設定した面の境界条件は

時間ステップが変わっても継続されているものなのでしょうか？

2014年9月11日木曜日 15時04分54秒 UTC+9 ohbuchi:

[ohbuchi]

各TimeディレクトリのTファイルを見ればわかる様に、各パッチの境界条件はそのまま引き継がれて
記述されますので、当然継続されます。
私の方で計算したところ、キャプチャ９に近い結果となっています。
ケースファイルをtarアーカイブで添付します。大きな違いは乱流計算になっている点です。
他は大きく変わったところはないはずです。
２つのロール状の循環流が、安定せず搖動するため定常解が得られていません。
尚、熱物性値を空気にすると対称な２つのロールになって1000ステップ程で収束します。
以上、ご参考まで。


2014年9月12日金曜日 7時33分02秒 UTC+9 hiro:

[hiro]

OpenFOAM学生のhiroです．

ohbuchi様，いつもありがとうございます．

ohbuchi様の添付ファイルthermCavity.tgzを参考に

実際の現象に近いものができました．

助かりました．

次のステップとして

ohbuchi様の添付ファイルで作成したものを

上面・下面の温度設定を逆にして計算を行い，

結果は対流が生じてしまいました．（キャプチャ10）

定常状態にもなりません（Time＝100000まで計算）

本来であれば，熱伝導のみで対流が起こらず，虹のようなきれいなものになるはずでした．

結果をみると

モデル枠（境界）の速度ベクトルでは常にX軸方向に向いています．

よって境界条件が何かしら悪影響しているのではと思います．

しかし，自分がみたところ

境界条件には問題が見られません．

そこで境界条件を記載せしますので，ご指摘頂けないでしょうか？

毎度助言を頂き，誠に恐縮ですが，よろしくお願いします．

------blockMeshDict-------

convertToMeters 0.001;
vertices
(
    (0  0  0)
    (80 0  0)

    (80 50 0)
    (0  50 0)
    (0  0  1)
    (80 0  1)
    (80 50 1)
    (0  50 1)

);
edges
(
);
blocks
(

    hex (0 1 2 3 4 5 6 7) (40 20 1) simpleGrading (1 1 1)

);
boundary
(
    fixedWalls
    {
        type wall;
        faces
        (

            (0 4 7 3)
            (1 2 6 5)

        );
    }

    cold
    {
        type wall;
        faces
        (
            (0 1 5 4)
        );
    }

    hot
    {
        type wall;
        faces
        (

            (3 7 6 2)
        );
    }
    empty
    {
        type  empty;
        faces
        (
            (0 3 2 1)

            (4 5 6 7)
        );
    }
);
mergePatchPairs
(
);

------0/T------

dimensions      [0 0 0 1 0 0 0];

internalField   uniform 293;

boundaryField
{
    fixedWalls
    {
        type            zeroGradient;
    }

    hot
    {
        type            fixedValue;
        value           uniform 313.15; // 40 degC
    }
    cold
    {
        type            fixedValue;
        value           uniform 288.15; // 15 degC
    }

    empty
    {
        type  empty;
    }
}

------0/U------

dimensions      [0 1 -1 0 0 0 0];

internalField   uniform (0 0 0);

boundaryField
{
    fixedWalls
    {
        type            fixedValue;
        value           uniform (0 0 0);
    }

    hot
    {
        type            fixedValue;
        value           uniform (0 0 0);
    }

    cold
    {
        type            fixedValue;
        value           uniform (0 0 0);
    }

    empty
    {
        type  empty;
    }
}

------0/alphat------

dimensions      [0 2 -1 0 0 0 0];

internalField   uniform 0;

boundaryField
{
    fixedWalls
    {
        type            alphatJayatillekeWallFunction;
        Prt             0.85;
        value            $internalField;
    }
    hot
    {
        type            alphatJayatillekeWallFunction;
        Prt             0.85;
        value           $internalField;
    }
    cold
    {
        type            alphatJayatillekeWallFunction;
        Prt             0.85;
        value           $internalField;
    }

    empty
    {
        type  empty;
    }
}

------0/epsilon------

dimensions      [0 2 -3 0 0 0 0];

internalField   uniform 0.01;

boundaryField
{
    fixedWalls
    {
        type            epsilonWallFunction;
        value            $internalField;
    }
    hot
    {
        type            epsilonWallFunction;
        value            $internalField;
    }
    cold
    {
        type            epsilonWallFunction;
        value            $internalField;
    }

    empty
    {
        type  empty;
    }
}

------0/k------

dimensions      [0 2 -2 0 0 0 0];

internalField   uniform 0.1;

boundaryField
{
    fixedWalls
    {
        type            kqRWallFunction;
        value            $internalField;
    }
    hot
    {
        type            kqRWallFunction;
        value            $internalField;
    }
    cold
    {
        type            kqRWallFunction;
        value            $internalField;
    }

    empty
    {
        type  empty;
    }
}

------0/nut------

dimensions      [0 2 -1 0 0 0 0];

internalField   uniform 0;

boundaryField
{
    fixedWalls
    {
        type            nutkWallFunction;
        value           uniform 0;
    }
    hot
    {
        type            nutkWallFunction;
        value           uniform 0;
    }
    cold
    {
        type            nutkWallFunction;
        value           uniform 0;
    }

    empty
    {
        type  empty;
    }
}

------0/p------

dimensions      [0 2 -2 0 0 0 0];

internalField   uniform 0;

boundaryField
{
    fixedWalls
    {
        type            calculated;
        value           $internalField;
    }

    hot
    {
        type            calculated;
        value           $internalField;
    }

    cold
    {
        type            calculated;
        value           $internalField;
    }

    empty
    {
        type  empty;
    }
}

------0/p_rgh------

dimensions      [0 2 -2 0 0 0 0];

internalField   uniform 0;

boundaryField
{
    fixedWalls
    {
        type            fixedFluxPressure;
        rho             rhok;
        value           $internalField;
    }

    hot
    {
        type            fixedFluxPressure;
        rho             rhok;
        value           $internalField;
    }

    cold
    {
        type            fixedFluxPressure;
        rho             rhok;
        value           $internalField;
    }

    empty
    {
        type  empty;
    }
}

2014年9月12日金曜日 11時16分01秒 UTC+9 ohbuchi:

[ohbuchi]

コンタのUのカラースケールを見て下さい。非常に小さい値で、ほぼ対流ゼロと

言ってよい結果です。

2014年9月18日木曜日 16時10分42秒 UTC+9 hiro:

[hiro]

ohbuchi様，ありがとうございます．

速度の数値がかなり小さいものになっていることに

気づきませんでした．

気づかせて頂き，ありがとうございます．

しかし，結果が定常にならない点や温度の結果が気になります．（キャプチャ12）

予想ではキャプチャ11のようになると思っています．

以上のことは境界条件，もしくは計算モデルが原因なのでしょうか？

2014年9月18日木曜日 21時35分47秒 UTC+9 ohbuchi:

[hiro]

すみません，添付ファイルをし忘れました．

2014年9月19日金曜日 15時58分03秒 UTC+9 hiro:

[ohbuchi]

おかしいですね。私の計算結果は安定成層した温度分布になっています。

添付画像を御覧ください。

最新のOpenFOAM-2.3.xで計算しました。計算ステップ数は２００００です。

お使いのバージョンにバグがあるのかも知れませんね。

2014年9月19日金曜日 16時00分05秒 UTC+9 hiro:

[hiro]

ohbuchi様

学生のhiroです．

先ほど，いろいろ試した結果

ohbuchi様の計算結果と同じ結果になりました．（キャプチャ1）

定常状態にならないのは変わりませんが...

先日の結果（キャプチャ12）になった原因は

fvSolutionの緩和係数Tの値を1にしていたことでした．

今回は緩和係数Tの値を0.5にした結果です．

以上から

自分のやりたい数値シミュレーションには

緩和係数Tが大きく影響することがわかりました．

どうもありがとうございます．

2014年9月20日土曜日 8時25分50秒 UTC+9 ohbuchi:

[ohbuchi]

バグでなく安心しました。

何をもって定常状態とみなすかだと思います。

安定成層温度分布が得られ、十分に小さいレベルの速度が残った状態であれば

完全に収束しなくとも、準定常状態であり、それを解とみなしても実用上の

問題は生じないことが多いでしょう。

2014年9月20日土曜日 20時37分05秒 UTC+9 hiro:

[hiro]

こんばんわ，学生のhiroです．

今回は自然対流シミュレーションにおける熱流束の設定についての

質問をさせて頂きます． 

今，二次元の自然対流シミュレーション（水）を行っており

上面冷却・底面加熱の境界条件に熱流束の設定を行いました．

計算も問題なく終えることができましたが，

なぜか結果における温度場では

温度差が1℃もない状態になっています．

（キャプチャ5,   プルームの形は理想的なのですが…）

今までは温度差が約10℃くらいある状態でした．

熱流束の設定では

上面では -50W/m2（15℃），底面では50W/m2（40℃）を想定しています．

あるスレッドには

「ついついqに与えたい熱流をそのまま書いてしまいがちですが、

非圧縮性ソルバを用いる場合は密度で割らなければならないことに注意が必要」

ということだったので水の密度で割り，それぞれ-50W/m2，50W/m2で

計算を行いましたが，温度場はどこかの国旗のようなになりました．（キャプチャ6）

熱流束の設定をするだけで温度差がなくなるものなのでしょうか？

それとも熱流束の設定に不備があったのでしょうか？

境界条件のほうを記載いたしますので

どなたかアドバイス等お願い致します．

------------0/T--------------

dimensions      [0 0 0 1 0 0 0];

internalField   uniform 288.15;

boundaryField
{
    fixedWalls
    {
        type            zeroGradient;
    }

    wall
    {
        type  empty;
    }

    hot
    {
        type            turbulentHeatFluxTemperature;
        heatSource      flux;        // flux (熱流束) または power (熱量)
        q               uniform 50; // 熱流束 [W/m2] または 熱量 [W]
        alphaEff        alphaEff;   // alphaEff field name;

        value           uniform 313.15; // 40 degC
    }

    cold
    {
        type            turbulentHeatFluxTemperature;
        heatSource      flux;        // flux (熱流束) または power (熱量)
        q               uniform -50; // 熱流束 [W/m2] または 熱量 [W]
        alphaEff        alphaEff;    // alphaEff field name;

        value           uniform 288.15;// 15 degC
    }

}

--------------U-------------------

dimensions      [0 1 -1 0 0 0 0];

internalField   uniform (0 0 0);

boundaryField
{
    fixedWalls
    {
        type            fixedValue;
        value           uniform (0 0 0);
    }

    wall
    {
        type  empty;
    }

    hot
    {
        type            fixedValue;
        value           uniform (0 0 0);
    }

    cold
    {
        type            fixedValue;
        value           uniform (0 0 0);
    }
}

--------------alphat-------------------

dimensions      [0 2 -1 0 0 0 0];

internalField   uniform 0;

boundaryField
{
    fixedWalls
    {
        type            alphatJayatillekeWallFunction;
        Prt             0.85;

        value           uniform 0;
    }

    wall
    {
        type  empty;
    }

    hot
    {
        type            alphatJayatillekeWallFunction;
        Prt             0.85;

        value           uniform 0;
    }

    cold
    {
        type            alphatJayatillekeWallFunction;
        Prt             0.85;

        value           uniform 0;
    }
}

--------------RASProperties-------------------

RASModel        laminar;

turbulence      on;

printCoeffs     on;

--------------nut-------------------

dimensions      [0 2 -1 0 0 0 0];

internalField   uniform 0;

boundaryField
{
    fixedWalls
    {
        type            nutkWallFunction;
        value           uniform 0;
    }

    wall
    {
        type  empty;
    }

    hot
    {
        type            nutkWallFunction;
        value           uniform 0;
    }

    cold
    {
        type            nutkWallFunction;
        value           uniform 0;
    }
}

--------------p-------------------

dimensions      [0 2 -2 0 0 0 0];

internalField   uniform 0;

boundaryField
{
    fixedWalls
    {
        type            calculated;
        value           $internalField;
    }

    wall
    {
        type  empty;
    }

    hot
    {
        type            calculated;
        value           $internalField;
    }

    cold
    {
        type            calculated;
        value           $internalField;
    }
}

何卒よろしくお願いします．

[ohbuchi]

温度の計算がうまくいっていない様です。

恐らくalphaEffが０になっており、解が初期値から変化していないと思われます。

層流で計算していることでalphaEffが正しく設定されていないのかも知れません。

laminar.HではalphaEffは正しくセットされている様に見えますが。

通常の勾配条件かgroovyBCを使って試すと良いでしょう。

2014年9月29日月曜日 18時28分51秒 UTC+9 hiro:

[hiro]

ohbuchi様のアドバイスを参考に

 turbulentHeatFluxTemperatureで境界条件を設定し，

0/TにalphaEffの値と通常の勾配条件（gradT）を書きこんで再度計算しました．

しかしコンパイルが通らず，計算にエラーが出て困っております．

ohbuchi様のアドバイスを間違って捉えてしまってるのでしょうか…

もしかして0/alphatを変更するのでしょうか？

0/Tの境界条件およびエラー内容を記載しますので

どなたかアドバイスお願いします．

------------0/T--------------

dimensions      [0 0 0 1 0 0 0];

internalField   uniform 288.15;

boundaryField
{
    fixedWalls
    {
        type            zeroGradient;
    }

    wall
    {
        type  empty;
    }

    hot
    {
        type            turbulentHeatFluxTemperature;
        heatSource      flux;        // flux (熱流束) または power (熱量)
        q               uniform 50; // 熱流束 [W/m2] または 熱量 [W]

        alphaEff        uniform 110e-6;   // alphaEff field name;
　　　 gradT=q/(rho0*Cp0)/alphaEff;

        value           uniform 313.15; // 40 degC
    }

    cold
    {
        type            turbulentHeatFluxTemperature;
        heatSource      flux;        // flux (熱流束) または power (熱量)
        q               uniform -50; // 熱流束 [W/m2] または 熱量 [W]

        alphaEff        uniform 110e-6;    // alphaEff field name;
 　　　gradT=q/(rho0*Cp0)/alphaEff;

        value           uniform 288.15;// 15 degC
    }

}

--------------エラー内容-----------------

Starting time loop

Time = 1

DILUPBiCG:  Solving for Ux, Initial residual = 1.39746e-15, Final residual = 1.39746e-15, No Iterations 0
DILUPBiCG:  Solving for Uy, Initial residual = 2.36945e-15, Final residual = 2.36945e-15, No Iterations 0

--> FOAM FATAL ERROR:

    request for volScalarField uniform from objectRegistry region0 failed
    available objects of type volScalarField are

9
(
rhok
alphaEff
p_rgh
nu
gh
alphat
p
T
p_rghPrevIter
)

    From function objectRegistry::lookupObject<Type>(const word&) const
    in file /home/opencfd/OpenFOAM/OpenFOAM-2.3.0/src/OpenFOAM/lnInclude/objectRegistryTemplates.C at line 198.

FOAM aborting

#0  Foam::error::printStack(Foam::Ostream&) in "/opt/openfoam230/platforms/linuxGccDPOpt/lib/libOpenFOAM.so"

#1  Foam::error::abort() in "/opt/openfoam230/platforms/linuxGccDPOpt/lib/libOpenFOAM.so"
#2  Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> const& Foam::objectRegistry::lookupObject<Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> >(Foam::word const&) const in "/opt/openfoam230/platforms/linuxGccDPOpt/lib/libfiniteVolume.so"
#3  Foam::incompressible::turbulentHeatFluxTemperatureFvPatchScalarField::updateCoeffs() in "/opt/openfoam230/platforms/linuxGccDPOpt/lib/libincompressibleTurbulenceModel.so"
#4  Foam::fvMatrix<double>::fvMatrix(Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> const&, Foam::dimensionSet const&) in "/opt/openfoam230/platforms/linuxGccDPOpt/bin/buoyantBoussinesqSimpleFoam"
#5 
 in "/opt/openfoam230/platforms/linuxGccDPOpt/bin/buoyantBoussinesqSimpleFoam"
#6  __libc_start_main in "/lib/i386-linux-gnu/libc.so.6"
#7 
 in "/opt/openfoam230/platforms/linuxGccDPOpt/bin/buoyantBoussinesqSimpleFoam"
中止 (コアダンプ)

どうぞよろしくお願いします．

2014年9月30日火曜日 6時00分49秒 UTC+9 ohbuchi:

[ohbuchi]

これはインラインコードのつもりでしょうか？
これでは構文が間違っているのでコンパイルは通りません。

私が述べたのは、普通のgroovyBCや通常の勾配型境界条件です。
groovyBCで過去ログを調べてみて下さい。

2014年10月3日金曜日 14時37分47秒 UTC+9 hiro:

[hiro]

こんばんわ，hiroです．

ohbuchi様が言っていることが

ようやく分かりました．

通常の勾配型境界条件というのは

fixedGradientのことを言っていたのですね．

このfixedGradientでフーリエの法則（q=gradT/k）から

温度勾配を計算して与えました．

ご教授，ありがとうございました．

現在，その境界条件（fixedGradient）で計算を行っています．

しかし，その境界条件について疑問があります．

今，2次元モデルで水に満ちている想定でシミュレーションを行っているのですが

その境界条件（fixedGradient）では，どの場所を境界と見ているのでしょうか？

たとえば

①メッシュの縁（キャプチャ10）

②メッシュのブロック（キャプチャ11）

③その他　　　　　　　　　　　　　　　　　　　どれなのでしょう？

ご存じであれば，教えてもらえないでしょうか．

よろしくお願いします．

2014年10月6日月曜日 12時41分56秒 UTC+9 ohbuchi:

[ohbuchi]

境界パッチはセルフェースで定義されますので、当然①になります。

2014年10月14日火曜日 22時59分20秒 UTC+9 hiro:

[hiro]

こんにちわ，hiroです．

ohbuchi様，いつもありがとうございます．

境界条件がどこで設定しているかでいろいろと傾向が

変わってくるので，ご教授頂いて助かりました．

今後とも，よろしくお願いします． 

[hiro]

こんばんわ，hiroです．

今回，流体密度の変更について質問させて頂きます．

詳しく言いますと，

今，水で行っている数値シミュレーションを

他の流体に変えて行えないかと考えております．

他の流体というのは，

ある温度で急激な密度変化が生じるものです．

この流体で自然対流シミュレーションを行うにあたり，

もしかするとbuoyantBoussinesqSimpleFoamソルバーの改造と

ソルバー内のcreateFields.H等を改造しなければならないのでしょうか？

まだソルバーの改造については勉強不足であるので

ご教授頂けたら幸いです．

どうぞ，よろしくお願いします．

[ohbuchi]

こんにちは。
ブシネ近似で扱えない程急激な密度変化の場合、buoyantSimpleFoamなどにすべき
でしょう。
蒸発を伴う様な場合は、simpleReactingParcelFoamなどが必要になると思います。
熱物性モデルで表現できる物性の場合、特にプログラムを修正する必要はありません。

2014年10月28日火曜日 22時35分26秒 UTC+9 hiro:

[hiro]

hiroです．

ohbuchi様，いつもお世話になっております．

 

急激な密度変化はbuoyantSimpleFoamということですが，

非圧縮性流体も取り扱う事が可能なのでしょうか？

ちなみに急激な密度変化といっても

20℃で密度：970kg/m3の流体が30℃くらいになると密度：950kg/m3と

約20kg/m3ほど低下するものです．

これがブジネスク近似で扱えないほど急激な変化とみるか

自分の解釈では計りかねております．

実際にブジネスク近似がどの程度まで許容できるのか

まだ理解できておりません．

（例えば，密度差・温度差がどの程度までいいのか）

→他のスレでは数十℃と明記されていましたが…

ご存知であれば，教えて頂けませんか？

よろしくお願いします．

[ohbuchi]

ブシネスク近似が成立するのはβ(T-T0)<<１が成立する範囲とされています。

単位温度変化当たりの密度変化が激しいとβが大きくなりますので要注意です。

またβが余り大きくなくとも、温度変化が大きければ上記成立条件を満たさなくなります。

http://ja.wikipedia.org/wiki/%E3%83%96%E3%82%B7%E3%83%8D%E3%82%B9%E3%82%AF%E8%BF%91%E4%BC%BC

2014年10月29日水曜日 22時44分25秒 UTC+9 hiro:

[hiro]

学生のhiroです．

ohbuchi様，貴重なアドバイスありがとうございます．

ブジネスク近似について，そのページはチェックをしていましたが

その数式には気づきませんでした．勉強になります．

あと先日にも申したのですが，

ある温度で急激な密度変化が生じるもので

自然対流シミュレーションを行いたいと思っており，

その流体密度を2つの温度の関数で表現したいと考えております．

例えば

30℃以下では関数①

30℃より高ければ関数②　と言ったものを

buoyantBoussinesqSimpleFoamのtransportPropertiesで

表示できないでしょうか？

[ohbuchi]

そういうことですと、ブシネスク近似は不適切だと思います。
浮力をρ0β（T-T0)と近似しているわけですから。
buoyantSimpleFoamにして、transportModelを何とかする必要があります。
標準で組み込まれている状態方程式では、多項式lで表現するicoPolynomialが
ありますが、それで表現できない場合は、カスタムモデルを作成するしかなさそうです。


2014年10月30日木曜日 17時38分58秒 UTC+9 hiro:

[hiro]

 

先日のohbuchi様のアドバイスから

buoyantSimpleFoamで水の自然対流シミュレーションを行っております．

（圧縮性ソルバで非圧縮性流体のシミュレーション）

そしてconstant/thermophysicalProperties で

水の物性値を用いて行っていますが，計算途中でエラーが出てしまいます．

おそらくthermophysicalPropertiesでの

多項式（polynomial，hPolynomial，icoPolynomial）の設定に

問題があると思っております．

そこで今回は物性値を多項式で設定することについて質問させて頂きます．

PENGUINITISのHPを参考に

水の物性値（密度，粘度，熱伝導率，比熱）を多項式で設定しております．

ここで設定しているのは，各物性値の多項式で用いられる係数を指定していることは分かります．

ただ，その物性値の多項式がどのような式なのか分かりません．

（いろいろと調べてはみたのですが…）

各物性値の多項式について

ご存知であれば教えて頂けないでしょうか？

それから以下に使用しているthermophysicalPropertiesを記載しますので

何か気になるところ・不備があればご指摘お願いします．

--------------------------constant/thermophysicalProperties -------------------------------

thermoType
{
    type            heRhoThermo;
    mixture         pureMixture;

    transport       polynomial;
    thermo          hPolynomial;
    equationOfState icoPolynomial;
    specie          specie;
    energy          sensibleEnergy;
}

mixture
{
    specie
    {
        nMoles          1;// モル数
        molWeight       18; // モル質量
    }
    equationOfState
    {
        rhoCoeffs<8>    ( 4.0097 -0.016954 3.3057e-05 -3.0042e-08 1.0286e-11 0 0 0 );//密度係数
    }

    thermodynamics
    {
        Hf              0; // 標準生成エンタルピー [J/kg]
        Sf              0; // 標準エントロピー [J/kg-K]
        CpCoeffs<8>     ( 948.76 0.39171 -0.00095999 1.393e-06 -6.2029e-10 0 0 0 );// 比熱係数
    }

    transport
    {
          muCoeffs<8>     ( 1.5061e-06 6.16e-08 -1.819e-11 0 0 0 0 0 ); //粘性係数係数
          kappaCoeffs<8>  ( 0.0025219 8.506e-05 -1.312e-08 0 0 0 0 0 ); //熱伝導率係数
    }
}

どうぞよろしくお願いします．

[hiro]

 一応，エラー内容も記載いたします．

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Create time

Create mesh for time = 0

Reading g
Reading thermophysical properties

Selecting thermodynamics package

{
    type            heRhoThermo;
    mixture         pureMixture;
    transport       polynomial;
    thermo          hPolynomial;
    equationOfState icoPolynomial;
    specie          specie;

    energy          sensibleEnthalpy;
}

Reading field U

Reading/calculating face flux field phi

Creating turbulence model

Selecting RAS turbulence model laminar
Calculating field g.h

Reading field p_rgh

No finite volume options present

Radiation model not active: radiationProperties not found
Selecting radiationModel none

SIMPLE: convergence criteria
    field p_rgh  tolerance 0.0001
    field U  tolerance 0.0001
    field h  tolerance 0.0001
    field "(k|epsilon|omega)"  tolerance 0.001

Starting time loop

Time = 1

DILUPBiCG:  Solving for Ux, Initial residual = 0.999998, Final residual = 0.00857706, No Iterations 1
DILUPBiCG:  Solving for Uy, Initial residual = 1, Final residual = 0.00316971, No Iterations 1
DILUPBiCG:  Solving for h, Initial residual = 1, Final residual = 0.00747037, No Iterations 1
GAMG:  Solving for p_rgh, Initial residual = 0.981306, Final residual = 0.00289608, No Iterations 5
time step continuity errors : sum local = 606775, global = -1.08388e-11, cumulative = -1.08388e-11

#0  Foam::error::printStack(Foam::Ostream&) in "/opt/openfoam230/platforms/linuxGccDPOpt/lib/libOpenFOAM.so"

#1  Foam::sigFpe::sigHandler(int) in "/opt/openfoam230/platforms/linuxGccDPOpt/lib/libOpenFOAM.so"
#2  Uninterpreted:
#3 
 in "/opt/openfoam230/platforms/linuxGccDPOpt/bin/buoyantSimpleFoam"
#4  __libc_start_main in "/lib/i386-linux-gnu/libc.so.6"
#5 
 in "/opt/openfoam230/platforms/linuxGccDPOpt/bin/buoyantSimpleFoam"
浮動小数点例外 (コアダンプ)

[ohbuchi]

polynomialTransportの定義は、

　f(T)=C0+C1*T+C2*T^2+...+C7*T^7

です。＜８＞で定義すると温度Tの７次式になります。

2014年11月7日金曜日 16時34分54秒 UTC+9 hiro:

[ohbuchi]

初回の連続の式残差がとても大きく、直ぐに発散していることから

熱物性の定義に無理がある様に見えますが、これだけでは何とも言えません。

緩和係数を下げて様子をみるのも良いかも知れません。

2014年11月8日土曜日 23時53分38秒 UTC+9 hiro:

[hiro]

ohbuchi様，いつも貴重なアドバイスありがとうございます．

ohbuchi様が熱物性の定義に無理があると言ったのは

多項式で用いられる係数の値についてでしょうか？

ここで記載してある係数はPENGUINITISで

記載してあるものをそのまま引用させて頂きました．

polynomialTransportの定義が

　f(T)=C0+C1*T+C2*T^2+...+C7*T^7　ということなので

再度，これに合わせて水の物性値に書き換えました．

しかし，計算は前と同じところで発散してしまい，

緩和係数の値を低くしてもエラーが出てしまいます．

この原因を考えたのですが，

多項式の定義で使用されるT（温度）に起因しているのではないでしょうか？

多項式で使用されているTが

絶対温度[℃]か摂氏温度[K]かで

物性値の値もかなり変動してくると思いました．

多項式のTは絶対温度・摂氏温度どちらなのでしょうか？

一応，（係数を）水の物性値に書き換えたthermophysicalPropertiesも

記載しますのでアドバイスよろしくお願いします．

--------------------------constant/thermophysicalProperties -------------------------------

thermoType

{
    type            heRhoThermo;
    mixture         pureMixture;
    transport       polynomial;
    thermo          hPolynomial;
    equationOfState icoPolynomial;
    specie          specie;
    energy          sensibleEnthalpy;
}

mixture
{
    specie
    {
        nMoles          1;
        molWeight       18;
    }

    equationOfState
    {
      rhoCoeffs<8> (999.1 0 0 0 0 0 0 0);
    }

    thermodynamics
    {
        Hf              0;
        Sf              0;
      CpCoeffs<8> (4187 0 0 0 0 0 0 0);

    }

    transport
    {
     muCoeffs<8> (1.154e-03 0 0 0 0 0 0 0);
     kappaCoeffs<8> (0.588 0 0 0 0 0 0 0);

    }
}

2014年11月9日日曜日 20時21分58秒 UTC+9 ohbuchi:

[ohbuchi]

圧縮性ソルバ（Boussinesqでないもの）では［K］です。

2014年11月11日火曜日 1時39分43秒 UTC+9 hiro:

[hiro]

hiroです．いつもお世話になっております．

前回の質問，ohbuchi様のご回答からいろいろと試してみました．

ohbuchi様，いつもアドバイスありがとうございます．

いろいろと試した結果，

やはりbuoyantBoussinesqSimpleFoamのソルバーを改造して

急激な密度変化を再現したいということに至りました．

そこで今回は，ソルバーの改造について質問させて頂きます．

改造したいソルバーは，buoyantBoussinesqSimpleFoam （定常熱流体解析ソルバ）．

上記ソルバーのcreateFields.Hを改造したと思っております．

具体的には，createFields.H の浮力項について （下記に記します）

--------------------------------------createFields.H--------------------------------------------

（省略）

    // Kinematic density for buoyancy force
    volScalarField rhok
    (
        IOobject
        (
            "rhok",
            runTime.timeName(),
            mesh
        ),
        1.0 - beta*(T - TRef)
    );

（省略）

上記の赤字部分を変更させて，急激な密度変化を再現したいと思います．

さらに言うと，赤字部分に式を2つ用いて if 文で表現ができないかと考えておりますが，

実際には，この自分の考えは可能なのでしょうか？

以上が質問の内容になります．

どのような意見でも構いませんので，よろしくお願いします．

2014年11月11日火曜日 4時43分20秒 UTC+9 ohbuchi:

[hiro]

学生のhiroです．

今回はクーラン数について お聞きしたいことがあります．

ユーザーガイドに記載されてあるクーラン数の式で，

|U|はセルを通る流速の大きさとありますが，

これはどのように決まるのでしょうか？

時間刻みΔtを決めるとき，この|U|が未知数だと計算ができませんよね？

自分で予想しろということでしょうか？

自分のやろうとシミュレーションは

閉鎖空間内の自然対流なので，それほど速い流速ではないと思いますが…

ご回答よろしくお願いします．

[ohbuchi]

まず、新しい質問は新しいスレッドを立ててお願いします。
あとでログ検索をする際に便利なためです。

時間刻みを自動調節する機能があります。
下記の様にcontrolDictに記述すると、最大Co数を満たす様にDtを自動調節します。

adjustTimeStep     on;
maxCo                  0.2;

ご参考まで。

2014年12月23日火曜日 3時51分19秒 UTC+9 hiro:

[hiro]

ohbuchi様，ご回答ありがとうございます．

早速やってみたいと思います．

それから

今後，質問する際は新しいスレッドを立てて行いたいと思います．

2014年12月23日火曜日 15時08分58秒 UTC+9 ohbuchi: