chtMultiRegionFoam解析について

[art.sugawara]

いつもお世話になっております。

熱解析の素人なのですが無謀にもchtMultiRegionFoamの解析にトライしています。
アルミ板上部に300K冷水を流し、下部は500K熱水を流します。アルミ板に並行した流れ
ではなくL字の流路にしています。設定などはchtMultiRegionFoamのチュートリアル、
multiRegionLiquidHeaterやsnappyMultiRegionHeaterなどを参考にしました。

OpenFOAM 2.2.2での解析結果が、OF222結果.pngです。同じ解析を商用CAE
（使っている人なら何かは一目瞭然ですが）で解析したのが商用結果.pngです。

両者の結果を見比べると流速ベクトルは最大値こそ10%程度の差がありますが、
どちらもそれなりです。問題はアルミ板内部の温度分布です。冷水／熱水が突き
当たる箇所が最低／最高温度を示す商用結果は理に適っているように思えます。
OpenFOAMの結果はそれとはかなり違いがあります。OpenFOAMモデルを添付
しましたので、どこに問題があるのか、その原因などのアドバイスをお願いします。

[ohbuchi]

熱物性値の条件が違うのではないでしょうか？

上下の流体領域のコーナー部では淀み域が形成され、むしろ出口に近い領域で

流速が早くなっています。仮に熱流束がインターフェースで一定だとすると

流れ方向に徐々に温度が上昇・下降していくOpenFOAMの結果は自然な気が

します。

2014年9月14日日曜日 23時07分18秒 UTC+9 art.sugawara:

[art.sugawara]

ohbuchi さん、いつもお世話になっております。

 

商用CAE結果の方が理にかなっていると思ったのは、流れが強制的に曲げられる

ことで、乱流の熱伝達が卓越するだろうと、すると冷水や熱水がぶつかる箇所での

熱移動が促進されると考えたからです。

 

ふと疑問に思い、chtMultiRegionFoamのログを眺めたら、

 

> Selecting turbulence model type RASModel
> Selecting RAS turbulence model realizableKE
> --> Upgrading k to employ run-time selectable wall functions
>     Backup original k to k.old
>     Writing updated k
> --> Upgrading epsilon to employ run-time selectable wall functions
>     Backup original epsilon to epsilon.old
>     Writing updated epsilon
> --> Creating mut to employ run-time selectable wall functions
>     Writing new mut
> --> Creating alphat to employ run-time selectable wall functions
>     Writing new alphat
> realizableKECoeffs

との記載がありました。参考にしたチュートリアルが層流だったため、alphatや

mutは入れ忘れてしまいました。ただ、これらはcompressibleな壁関数の設定

なので適切か吟味が必要かもしれません。しばらく格闘してみます。

2014年9月17日水曜日 22時53分45秒 UTC+9 ohbuchi:

[art.sugawara]

いつもお世話になっております。

 

なかなか時間が取れずに検討途中のままでしたが、最近ちょっと
時間が作れたので再検討していますが、難航しています。

 

chtMultiRegionFoamのチュートリアルを参考にしてモデルを
作成、さらにソルバー実行時に出るメッセージから、alphatと
mutが不足していると考えて0フォルダに追加、systemフォルダ
のchangeDictionaryDictも修正しました。

 

添付の残差プロット（heat_ex2b#1_residuals.png）から分か
るように流れ解が収束していない段階の結果ですが、流速ベク
トル（heat_ex2b#1.png）がおかしいです。特に流入口付近。
時刻歴で見ると広がったりすぼまったり振動しています。

 

問題は
system/bottomWater/changeDictionaryDict
system/topWater/changeDictionaryDict
のalphat、mutの流入口（あるいは流出口も？）設定にあると
睨んでいます。もしかしたらepsilonやkも問題なのかも
しれません。

 

悲しいことにチュートリアルの設定と比較しても私のスキルでは
問題箇所が見いだせずにいます。モデルを添付しましたので
ぜひご指摘をお願いいたします。

 

system/bottomWater/changeDictionaryDictのみ以下に
コピペします。
----------------------------------
/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  2.2.2                                 |
|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    object      changeDictionaryDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

dictionaryReplacement
{
    T
    {
        internalField   uniform 500;

        boundaryField
        {
            ".*"
            {
                type            zeroGradient;
            }

            minZ
            {
                type            fixedValue;
                value           uniform 500;
            }

            maxX
            {
                type            inletOutlet;
                inletValue      uniform 500;
                value           uniform 500;
            }

            "bottomWater_to_.*"
            {
                type            compressible::turbulentTemperatureCoupledBaffleMixed;
                neighbourFieldName T;
                kappa           fluidThermo;
                kappaName       none;
                value           uniform 500;
            }
        }
    }

    U
    {
        internalField   uniform (0 0 0);

        boundaryField
        {
            ".*"
            {
                type            fixedValue;
                value           uniform (0 0 0);
            }

            minZ
            {
                type            fixedValue;
                value           uniform ( 0 0 0.01 );
            }

            maxX
            {
                type            inletOutlet;
                inletValue      uniform ( 0 0 0 );
                value           uniform ( 0.01 0 0 );
            }
        }
    }

    alphat
    {
        internalField   uniform 0;

        boundaryField
        {
            ".*"
            {
                type            compressible::alphatWallFunction;
                Prt             0.85;
                value           uniform 0;
            }

            minZ
            {
                type            fixedValue;
                value           uniform 0;
            }

            maxX
            {
                type            inletOutlet;
                inletValue      uniform 0;
                value           uniform 0;
            }
        }
    }

    epsilon
    {
        internalField   uniform 0.01;

        boundaryField
        {
            ".*"
            {
                type            compressible::epsilonWallFunction;
                value           uniform 0.01;
            }

            minZ
            {
                type            fixedValue;
                value           uniform 0.01;
            }

            maxX
            {
                type            inletOutlet;
                inletValue      uniform 0.01;
                value           uniform 0.01;
            }
        }
    }

    k
    {
        internalField   uniform 0.1;

        boundaryField
        {
            ".*"
            {
                type            compressible::kqRWallFunction;
                value           uniform 0.1;
            }

            minZ
            {
                type            fixedValue;
                value           uniform 0.1;
            }

            maxX
            {
                type            inletOutlet;
                inletValue      uniform 0.1;
                value           uniform 0.1;
            }
        }
    }

    mut
    {
        internalField   uniform 0;

        boundaryField
        {
            ".*"
            {
                type            mutkWallFunction;
                Cmu             0.09;
                kappa           0.41;
                E               9.8;
                value           uniform 0;
            }

            minZ
            {
                type            fixedValue;
                value           uniform 0;
            }

            maxX
            {
                type            inletOutlet;
                inletValue      uniform 0;
                value           uniform 0;
            }
        }
    }

    p_rgh
    {
        internalField   uniform 1e5;

        boundaryField
        {
            ".*"
            {
                type            fixedFluxPressure;
                value           uniform 1e5;
            }

            maxX
            {
                type            fixedValue;
                value           uniform 1e5;
            }
        }
    }

    p
    {
        internalField   uniform 1e5;

        boundaryField
        {
            ".*"
            {
                type            calculated;
                value           uniform 1e5;
            }

            maxX
            {
                type            calculated;
                value           uniform 1e5;
            }
        }
    }
}

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

 

2014年9月19日金曜日 18時04分35秒 UTC+9 art.sugawara:

[ohbuchi]

mutやalphatは内部で計算されるので、流出入境界ではtype calculatedで指定しましょう。


2014年10月3日金曜日 13時20分52秒 UTC+9 art.sugawara:

[art.sugawara]

ohbuchi 様、いつもありがとうございます。

 

alphatとmutの流出入口のtypeをcaluculatedに変えて流してみました。

 

流入口は流速指定なので0.1m/sという一定値になるはずですが、
添付グラフのようにbottomWater流入口のセルの流速が振動しています。
流速ベクトル図（ピンクに選択されているセルの時刻歴取得）も添付しました。

 

再度changeDictionaryDictを見直ししています。設定漏れを1点見つけたので

その結果が出たら再報告させていただきます。

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

[art.sugawara]

自己レスです。流入口の流速指定は0.1m/sではなく0.01m/sです。 

流速の振動は0.01m/sの約2倍から1/2倍の範囲で振動しています。

 

[ohbuchi]

blockMeshDictを拝見しましたが、maxX, minZがwallタイプになっているのが気になりました。リージョン分割後に
boundaryファイルは正しくpatchタイプに変更しましたか？


2014年10月8日水曜日 14時34分21秒 UTC+9 art.sugawara:

自己レスです。流入口の流速指定は0.1m/sではなく0.01m/sです。 

流速の振動は0.01m/sの約2倍から1/2倍の範囲で振動しています。

 

[art.sugawara]

ohbuchi 様、いつもお世話になっております。

見逃していた設定がご指摘の境界条件でした。changeDictionaryDictを修正して、
流出入口のタイプがpatchになるように修正してみましたが、流入口での流速が
振動する現象は解消されませんでした。引き続き、その他の流出入口境界条件
を種々試行中です。

2014年10月9日木曜日 15時29分12秒 UTC+9 ohbuchi:

[art.sugawara]

いつもお世話になっております。

 

境界条件を種々振ってみて現象を確認してみました。いままでの

投稿と重複しますが以下説明します。

典型的な流体のchangeDictionaryDictを末尾にコピペしました。

（分かりやすくするため流出入境界名をinlet, outletに書き換えています）

 

これで解析した結果のポスト図ファイルを添付します。（heat_ex2b#6.png）

問題は中央のアルミ板に下から熱水、上から冷水を吹き付けているのですが、

例えば熱水流れの下流の方がアルミ温度が高いという変な結果です。

 

changeDictionaryDictの流入口のepsilonとkをチュートリアル

multiRegionLiquidHeaterのようにfixedValue（以下参照）とすると、

熱水流れの上流のアルミ温度が高くなり理にかなっていると思いますが、

流速規定しているinletでの流速が時間的に振動する不具合発生します。

（heat_ex2b#7_no.8.png）

 

----------
epsilon
----------
inlet

{
    type            fixedValue;
    value           uniform 0.01;
}

----------
k
----------
inlet

{
    type            fixedValue;
    value           uniform 0.1;
}

----------

 

heat_ex2b#6の解析結果でbottomWaterのみの種々結果図を

添付いたしました。原因究明などに結びつく何かが見つかると

うれしいのですが皆様のアドバイスをお待ちしています。

 

changeDictionaryDict_#6

/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  2.2.2                                 |
|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    object      changeDictionaryDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

dictionaryReplacement
{
    boundary
    {
        inlet
        {
            type            patch;
        }
        outlet
        {
            type            patch;
        }
    }

    T
    {
        internalField   uniform 500;

        boundaryField
        {
            ".*"
            {
                type            zeroGradient;
            }

            inlet

            {
                type            fixedValue;
                value           uniform 500;
            }

            outlet
            {
                type            zeroGradient;
            }

            "bottomWater_to_.*"
            {
                type            compressible::turbulentTemperatureCoupledBaffleMixed;
                neighbourFieldName T;
                kappa           fluidThermo;
                kappaName       none;
                value           uniform 500;
            }
        }
    }

    U
    {
        internalField   uniform (0 0 0);

        boundaryField
        {
            ".*"
            {
                type            fixedValue;
                value           uniform (0 0 0);
            }

            inlet

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

            outlet
            {
                type            zeroGradient;
            }
        }
    }

    alphat
    {
        internalField   uniform 0;

        boundaryField
        {
            ".*"
            {
                type            compressible::alphatWallFunction;
                Prt             0.85;
                value           uniform 0;
            }

            inlet
            {
                type            calculated;
                value           uniform 0;
            }

            outlet
            {
                type            calculated;
                value           uniform 0;
            }
        }
    }

    epsilon
    {
        internalField   uniform 0.01;

        boundaryField
        {
            ".*"
            {
                type            compressible::epsilonWallFunction;
                value           uniform 0.01;
            }

            inlet
            {
                type            compressible::turbulentMixingLengthDissipationRateInlet;
                mixingLength    0.005;
                value           uniform 0.01;
            }

            outlet
            {
                type            zeroGradient;
            }
        }
    }

    k
    {
        internalField   uniform 0.1;

        boundaryField
        {
            ".*"
            {
                type            compressible::kqRWallFunction;
                value           uniform 0.1;
            }

            inlet
            {
                type            turbulentIntensityKineticEnergyInlet;
                intensity       0.05;
                value           uniform 0.1;
            }

            outlet
            {
                type            zeroGradient;
            }
        }
    }

    mut
    {
        internalField   uniform 0;

        boundaryField
        {
            ".*"
            {
                type            mutkWallFunction;
                Cmu             0.09;
                kappa           0.41;
                E               9.8;
                value           uniform 0;
            }

            inlet
            {
                type            calculated;
                value           uniform 0;
            }

            outlet
            {
                type            calculated;
                value           uniform 0;
            }
        }
    }

    p_rgh
    {
        internalField   uniform 1e5;

        boundaryField
        {
            ".*"
            {
                type            fixedFluxPressure;
                value           uniform 1e5;
            }

            outlet

            {
                type            fixedValue;
                value           uniform 1e5;
            }
        }
    }

    p
    {
        internalField   uniform 1e5;

        boundaryField
        {
            ".*"
            {
                type            calculated;
                value           uniform 1e5;
            }

            outlet

            {
                type            calculated;
                value           uniform 1e5;
            }
        }
    }
}

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

[nakagawa]

いくつか気になった点について，コメントさせて頂きます。

「流速規定しているinletでの流速が時間的に振動する不具合」
　これは，表示方法の問題だと推測します。ソルバの不具合ではなく。
　以前の投稿で，セルの速度をプロットされていました。今回，計算領域内部の速度ベクトルを見て，判断されていると思います。
　境界条件として既定しているのは，あくまでも流入面での速度です。patchでの速度は，既定値で一定だと思います。
　一致位置で指定した面に隣接するセルでの計算値が，振動しているのではないでしょうか。

　これは，振動が起こるような領域に境界面を設定していることに問題があるのではないでしょうか？
　計算領域を拡張して，速度を一定にしても問題ない位置まで，境界を遠ざける必要があるのでは。
　（商用ではOKなのに？という意見もあるでしょうが．．．）

　ベクトル表示方法も，少し気になります．cellCenterフィルタを使わずに，接点での速度ベクトルを表示しているのでしょうか。それだと，真の計算結果を表示できない可能性があると思います。とくに，メッシュが粗く，境界近くだと。


　また，メッシュが粗すぎることも問題の一因のように感じます。
　流路中央部，かなり粗いメッシュですよね。
　メッシュを少し細かくするとどうなるのでしょう。2次元計算でも良いのではと思います。
　比較されている商用ソフトのモデルとは，セル数は同じでしょうか？

　流速が0.01m/sとかなり遅いですが，乱流モデルを使って計算する必要があるのでしょうか。モデルの大きさが分からないので，何とも言えないところですが。
　層流で計算しても，同様な問題が発生するのでしょうか。

・結果表示位置について
　結果の図は，何の表示なのでしょうか？　y方向のどこかのスライス面ですか？それとも，計算領域を外から眺めているだけでしょうか？
　流体の温度がほぼ500Kとなっています。これはなぜなのでしょう？


ざっと，気になった点をコメントいたしました。雑ですみません。
既にご検討済みの点などございましたら，無視ししていただいて結構です。

なかがわ

2014年10月21日火曜日 14時58分41秒 UTC+9 art.sugawara:

[art.sugawara]

なかがわ 様

アドバイスいただきありがとうございます。
ご指摘内容を確認してみました。

> 「流速規定しているinletでの流速が時間的に振動する不具合」
> 　これは，表示方法の問題だと推測します。ソルバの不具合ではなく。
> 　以前の投稿で，セルの速度をプロットされていました。今回，計算領域内部の速度ベ
> クトルを見て，判断されていると思います。
> 　境界条件として既定しているのは，あくまでも流入面での速度です。patchでの速度
> は，既定値で一定だと思います。
> 　一致位置で指定した面に隣接するセルでの計算値が，振動しているのではないでしょ
> うか。

そういう観点が欠落しておりました。paraFoamでpatch面での流速時刻歴
グラフを作ったところ設定どおりの0.01m/sでした。
流入口patchでの他物理量も調べたところ、alphat、mut、p、p_rghに
振動が見られました。これがセル値流速のあばれ（振動）を引き起こして
いるのではないかと思われます。
もしこれが正しければ、

> 　これは，振動が起こるような領域に境界面を設定していることに問題があるのではな
> いでしょうか？
> 　計算領域を拡張して，速度を一定にしても問題ない位置まで，境界を遠ざける必要が
> あるのでは。
> 　（商用ではOKなのに？という意見もあるでしょうが．．．）

おっしゃるように流入口をもっと離してあげれば合理的な解になる可能性
があります。

流速の振動が起こらないケース（heat_ex2b#6）で流入口patchでの結果
を調べたら、epsilon=7.55e-9、k=3.75e-7と非常に小さくなっています。
もしかしたらこのせいで伝熱がうまくいかないのかもしれません。

つまり、
(1)振動はないが温度が？のモデル（heat_ex2b#6）では流入口の乱流設定を見直す。
(2)振動するが温度はそれなりのモデル（heat_ex2b#7_no.8）では流入口をもっと遠ざける。
という2通りが考えられるのかなと思っていますがいかがでしょうか？

> 　ベクトル表示方法も，少し気になります．cellCenterフィルタを使わずに，接点での
> 速度ベクトルを表示しているのでしょうか。それだと，真の計算結果を表示できない可
> 能性があると思います。とくに，メッシュが粗く，境界近くだと。

cellCenterは初耳です。セル中心／セル表面を意識していませんでした。

> 　また，メッシュが粗すぎることも問題の一因のように感じます。
> 　流路中央部，かなり粗いメッシュですよね。
> 　メッシュを少し細かくするとどうなるのでしょう。2次元計算でも良いのではと思います。
> 　比較されている商用ソフトのモデルとは，セル数は同じでしょうか？

商用ソフトではアルミ板がかなり粗くて全体で約17000ですが、OpenFOAMは
全領域同じメッシュ密度でセル数は約32000です。流体部のメッシュ密度は
同程度です。メッシュ数をもっと増やす検討はしてみます。

> 　流速が0.01m/sとかなり遅いですが，乱流モデルを使って計算する必要があるのでし
> ょうか。モデルの大きさが分からないので，何とも言えないところですが。
> 　層流で計算しても，同様な問題が発生するのでしょうか。

モデルはXとYの寸法が1m、Z寸法は1.25mです。管路の代表寸法が0.25m
なので、レイノルズ数を計算すると1000*0.01*0.25/0.0008=3125という
ことで乱流域と考えました。

> ・結果表示位置について
> 　結果の図は，何の表示なのでしょうか？　y方向のどこかのスライス面ですか？それ
> とも，計算領域を外から眺めているだけでしょうか？

y方向中央断面での表示です。

> 　流体の温度がほぼ500Kとなっています。これはなぜなのでしょう？

熱水はアルミ板近傍で若干の熱交換するのみで、ほとんどそのまま
流出してしまうのでそうなったのだと思います。

 

いただいたコメント・アドバイスを元にまた試行錯誤してみます。
今後もよろしくお願いいたします。

2014年10月21日火曜日 20時54分52秒 UTC+9 nakagawa:

[nakagawa]

商用ソフトではアルミ板がかなり粗くて全体で約17000ですが、OpenFOAMは
全領域同じメッシュ密度でセル数は約32000です。流体部のメッシュ密度は
同程度です。メッシュ数をもっと増やす検討はしてみます。

　そうですね，固体は粗くても良いかもしれませんね。　ただ，流体部はかなり粗く感じます。

　商用の解では，すみ部や流路壁面での逆流領域を再現しているようですが，今回のメッシュではとても無理そうです。

> 　流体の温度がほぼ500Kとなっています。これはなぜなのでしょう？
熱水はアルミ板近傍で若干の熱交換するのみで、ほとんどそのまま
流出してしまうのでそうなったのだと思います。

　固体壁面自体はもっと低温になっているようなので，水温も，もう少し低下しても良いのかと思いました。

　paraviewでの可視化で，現在はポイントでの温度表示を選択されていますが，この表示も，実際の計算結果を把握しにくくすると思います。特にセルが粗い場合など，表示上はなめらかな温度分布のように見せかけてしまいます。まずは，セル表示を選択されることをおすすめします。その表示でみると，壁面近傍第一セルでの温度も，ほとんど500Kとなるのでは，と危惧します。

おおぶち様のご指摘で、基本的な問題はクリアされた。後はメッシュが粗すぎる，という状態だとよいですね。メッシュ修正だけで良い結果がでることを祈っています。

壁関数を使うにしても，壁面近傍セルが大きすぎるのかも。ひとまず層流で試す方法もあるのではないでしょうか。

なかがわ 

2014年10月22日水曜日 15時04分32秒 UTC+9 art.sugawara:

[art.sugawara]

いつもお世話になっております。

 

すっきりと解決したわけではありませんが、すこし状況が整理
できて見えてきたものもありますので報告させていただきます。
添付図は流入口でのepsilon、kを種々振った結果です。

 

typeA_#9f.png：チュートリアルと同様にepsilon=0.01、k=0.1
流速規定の流入条件ですが流入直後に流れが乱れてしまいます。
流れ挙動はメッシュ密度を変えても流入導入部を長めにしても
変わりませんでした。

 

typeB_#8.png/typeB_#9d.png：以下のサイトを参考にして設定
http://www.geocities.jp/penguinitis2002/study/OpenFOAM/bc_settings.html
解析結果によれば流入口のepsilon=1e-9、k=4e-7程度でした。
流入口近傍での流速乱れは起きていません。出口流速は山なり
の分布なので層流に近いと思われます。
固体部の温度分布はメッシュ依存（#8は粗、#9dは精細）で
#8解は信頼置けないと思われます。

 

流入口のepsilonをmixingLengthで指定したtypeBはメッシュ密度
を高めることで信じるに足る解になったのではと思われます。

typeB_#9dがもっともらしいのかなというのが現時点の感触です。

2014年10月23日木曜日 2時14分11秒 UTC+9 nakagawa: