【OpenFOAM】球体周りの自然対流

物理・プログラミング・Python、Web制作・自己啓発

2023年11月23日 23:25

こんにちは。
球体周りの自然対流をやっていてなかなかうまくいかなかったのでメモとして残しておきます。
最終的にはこのような流れになります。

自然対流での問題

PENGUINITISさんの記事にも書かれているように自然対流で大気開放で流れが生じる場合にうまく流れが上に抜けていかないという問題があります。

こちらにも議論が上がっています。

まずチュートリアルをそのまま計算します。
上面が壁なので流れは壁で流速が0m/sなので旋回しているのが確認できます。

これを大気開放にすると逆流が発生するとのことで以下のように変更すると良いと・・・っでやってみた。

0/p_rgh

    ceiling
    {
//        type            fixedFluxPressure;
//        value           uniform 100000;
        type            prghPressure;
        p               uniform 100000;
        value           uniform 100000;
    }

あれ？
これでいいのか？
あと、乱流モデルなしにするとエラーになる。
OpenFOAMのチュートリアルは概ね計算さえできたらいいみたいなところがあるのでこの辺はメッシュの設定など適当かもしれない。

このチュートリアルについてはちゃんと検証する気はなかったので、圧力の境界条件に注意が必要ということを認識する程度にして、本当にやりたい球体周りの自然対流の計算に移ります。

球体周りの自然対流の設定

結局うまくいった設定を載せておきます。

メッシュはこのようになっています。
まず右の絵のようにblockMeshで球体周りのメッシュを作成して、上側の領域をextrudeMeshで広げました。

メッシュの作成方法はこちらにも記載しています。

15000ステップまで計算させてようやく流れが上に抜ける落ち着いた自然対流っぽい流れになりました。

0/U

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

internalField   uniform (0 0 0);

boundaryField
{
    Zmin
    {
        type            noSlip;
    }

    Zmax
    {
        type            pressureInletOutletVelocity;
        value           uniform (0 0 0);
    }

    "(Xmin|Xmax|Ymin|Ymax)"
    {
        type            symmetry;
    }

    ball
    {
        type            noSlip;
    }
}

0/p_rgh

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

internalField   uniform 100000;

boundaryField
{
    Zmin
    {
        type            fixedFluxPressure;
        value           uniform 100000;
    }

    Zmax
    {
//        type            fixedFluxPressure;
//        value           uniform 100000;
        type            prghPressure;
        p               uniform 100000;
        value           uniform 100000;
    }

    "(Xmin|Xmax|Ymin|Ymax)"
    {
        // type            fixedFluxPressure;
        // value           uniform 100000;
       type symmetry;
    }

    ball
    {
        type            fixedFluxPressure;
        value           uniform 100000;
    }
}

0/T
球体周りの温度は330Kとしました。

dimensions      [0 0 0 1 0 0 0];

internalField   uniform 293;

boundaryField
{
    Zmin
    {
        type            fixedValue;
        value           uniform 293.0;
    }

    Zmax
    {
//        type            fixedValue;
//        value           uniform 300.0;
        type            inletOutlet;
        inletValue      uniform 293;
        value           uniform 293;

    }

    "(Xmin|Xmax|Ymin|Ymax)"
    {
        // type            zeroGradient;
        type symmetry;
    }

    ball
    {
        type            fixedValue;
        value           uniform 330.0;
    }
}

ここからは乱流モデル用(標準k-ε)の設定
チュートリアルそのままなので高レイノルズ数型の壁関数などとしているが、これでいいのかは検証していない。

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

internalField   uniform 0.1;

boundaryField
{
    Zmin
    {
        type            kqRWallFunction;
        value           uniform 0.1;
    }

    Zmax
    {
        // type            kqRWallFunction;
        // value           uniform 0.1;
        type            inletOutlet;
        inletValue      $internalField;
        value           $internalField;
    }

    "(Xmin|Xmax|Ymin|Ymax)"
    {
        // type            kqRWallFunction;
        // value           uniform 0.1;
               type symmetry;

    }

    ball
    {
        type            kqRWallFunction;
        value           uniform 0.1;
    }
}

0/epsilon

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

internalField   uniform 0.01;

boundaryField
{
    Zmin
    {
        type            epsilonWallFunction;
        value           uniform 0.01;
    }

    Zmax
    {
        // type            epsilonWallFunction;
        // value           uniform 0.01;
        type                inletOutlet;
        inletValue          $internalField;
        value               $internalField;
    }

    "(Xmin|Xmax|Ymin|Ymax)"
    {
        // type            epsilonWallFunction;
        // value           uniform 0.01;
               type symmetry;

    }

    ball
    {
        type            epsilonWallFunction;
        value           uniform 0.01;
    }
}

0/nut（渦粘性）
kとεから計算される

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

internalField   uniform 0;

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

    Zmax
    {
        // type            nutkWallFunction;
        // value           uniform 0;
        type            calculated;
        value           uniform 0;
    }

    "(Xmin|Xmax|Ymin|Ymax)"
    {
        // type            nutkWallFunction;
        // value           uniform 0;
               type symmetry;

    }

    ball
    {
        type            nutkWallFunction;
        value           uniform 0;
    }
}

0/alphat（乱流拡散）

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

internalField   uniform 0;

boundaryField
{
    Zmin
    {
        type            compressible::alphatWallFunction;
        value           uniform 0;
    }

    Zmax
    {
        // type            compressible::alphatWallFunction;
        // value           uniform 0;
        type                calculated;
        value               $internalField;
    }

    "(Xmin|Xmax|Ymin|Ymax)"
    {
        // type            compressible::alphatWallFunction;
        // value           uniform 0;
               type symmetry;

    }

    ball
    {
        type            compressible::alphatWallFunction;
        value           uniform 0;
    }
}

熱伝達率の出力

熱伝達率の出力は以下で行えます。

heatTransferCoeff1
{
    // Mandatory entries (unmodifiable)
    type            heatTransferCoeff;
    libs            (fieldFunctionObjects);
    field           T;//<field>;
    patches         (ball);//(<patch1> <patch2> ... <patchN>);
    htcModel        fixedReferenceTemperature;//<htcModel>;
    //UInf            (0 0 0);
    //Cp              CpInf;
    //CpInf           1000;
    //rho             rhoInf;
    //rhoInf          1.2;
    TRef              293;

    // Optional (inherited) entries
    //result          <fieldResult>;
    region          region0;
    enabled         true;
    log             true;
    timeStart       0;
    timeEnd         20000;
    executeControl  timeStep;
    executeInterval 500;
    writeControl    timeStep;
    writeInterval   500;
}

球体周りの自然対流の文献値はありますので計算ツールで計算して比較してみましょう。

平均させると4.4[W/m2 K]。
大外れではなく、だいたいあっているような気がします。
しかし、球体なのに方向性がある局所熱伝達率のの分布になっているのが気になります。

熱伝達率はだいたい文献値通りになった。
しかし、こういう分布でいいのか？

↓計算ツールhttps://t.co/lAvAybzxdF pic.twitter.com/QhH3xa4F0T
— カマキリ🐲CAE頑張る (@t_kun_kamakiri) November 21, 2023

あと、今回は流速のinternalFieldを(0 0 0)でしましたが、球体周りはz方向に0.1m/sくらい出ているので、internalFieldを(0 0 0.1)で計算開始させた方が収束までが早いかもしれません。
定常解析なので初期状態に依存しないとはいえ、正解に近い状態から始めた方がいいですよね。

OpenFOAMをはじめるなら

OpenFOAMをはじめるなら手にしておきたい日本語の書籍がこちらです。
OpenFOAMをインストールしてこれから深く学んでいこうという方や、使い慣れてもレベルアップのために何回も読みたくなるくらい参考になる内容です。

初心者の方にはチュートリアルを動かしながら学んでいくのが良いでしょう。こちらはOpenFOAMをインストールしてチュートリアルを少し触っていくところからはじめて最後に重合格子をやるなど面白い内容です。

Twitter➡@t_kun_kamakiri
ブログ➡宇宙に入ったカマキリ（物理ブログ）
youtube➡理系ブログ

この記事が気に入ったらサポートをしてみませんか？