【OpenFOAM】熱流体のチュートリアルで発散した。改善の着眼点。
最近、OpenFOAMで熱流体をやりはじめました。
チュートリアルを適当に持ってきて色々計算させようと思ったけど、発散して計算が上手く進まなかったのでメモとして残しておきます。
熱流体、発散するから
— カマキリ🐲技術士勉強中 (@t_kun_kamakiri) January 19, 2023
重力無しにすると発散せず
乱流モデル入れると発散せず
ここから言えることは何だろ。
Twitterで適当に質問をすると意見をもらえることがあります。今回も適当に呟いたところ、気軽に返答されたので少し解決のヒントになりました。
OpenFOAMに関する質問でしたら以下のGoogle Groupsに投げかけると、有志で誰かが答えてくれますのでそちらの方がきっちり解決したい方には良いかもしれません。
ただし、質問をする場合は上手にしましょう。
・目的（何がやりたいか）
・解析条件（どこの設定に問題がありそうか）
手っ取り早いのは誰かが手元で計算実行させてエラー文を見て修正できるように、モデル一式を圧縮して添付することですかね。質問が丸投げで「あとはよろしく！」って人がいますが、無料でも有料関係なく状況をきっちり伝えないと答えようがありません。そのような丸投げはやめましょう。
Twitterは適当に呟くと適当に回答もらえるくらいの温度感なので、それでも良い方にはちょうど良い温度感でやり取りできて楽です。
ということで、熱流体をやってて「なんか、計算が進まないなー」って思ってTwitterで適当に呟いて、得られた回答をヒントに試した内容をメモとして残しておきます。
OpenFOAM-v2012
上手くいかないときは原因を考えて影響因子を絞る
計算が発散して「あれ？」って思った経緯は、チュートリアルのモデルが乱流モデル（RANS k-ε）を使っており、「乱流モデルの設定を考えるのが面倒なので、乱流モデル無しでしよう」と思って計算したら発散したからです。
自分が解析したい解析モデルの前にチュートリアルでテスト計算していたら、思わぬ発散が起こってしまったため色々と調べて分かったことを本記事でまとめておこうと思います。
発散する直前の計算ログ。
DILUPBiCGStab: Solving for Ux, Initial residual = 0.072134, Final residual = 0.00064676, No Iterations 1
DILUPBiCGStab: Solving for Uy, Initial residual = 0.0745978, Final residual = 0.000783793, No Iterations 1
DILUPBiCGStab: Solving for Uz, Initial residual = 0.136694, Final residual = 0.000918327, No Iterations 1
DILUPBiCGStab: Solving for h, Initial residual = 0.265034, Final residual = 0.00251917, No Iterations 1
DICPCG: Solving for p_rgh, Initial residual = 0.793862, Final residual = 0.00646339, No Iterations 20
time step continuity errors : sum local = 0.618086, global = 6.47959e-16, cumulative = 1.93099e-16
rho min/max : -5.40042 7.88036
ExecutionTime = 1.33 s ClockTime = 2 s
Time = 4continuity errors : sum local = 0.618086, この値が大きすぎることに気付きます。これは連続の式の誤差が大きすぎることを意味しており、圧力の残差もおそらく相当大きいでしょう。
完全に連続の式を満たさないにしても、感覚的にはcontinuity errors : sum local =1e-6~1e-7くらいのオーダーであれば正常です。場合によっては0.001くらいから徐々に小さくなってくる場合もありますが、今回のような単純なモデルで0.6というのは非常に大きな値であるという感じがします。
おそらく原因は「メッシュが粗いため局所的に流速変化が大きくなる部分で方程式を満たさず、誤差が溜まって発散」と考えられます。
これを仮定とおいて影響因子を割り出しスタディしました。
グラフホフ数($${\frac{gβΔTh^3}{\nu^2}}$$)を計算すると$${10^{10}}$$と自然対流で乱流になっていてもおかしくないので、乱流なのにメッシュが粗くて速度変化を再現できず誤差が溜まって発散したと思われます。
なので、乱流モデルを導入して乱流を数値的にモデル化したりするなど対応が必要です。今回は乱流モデルを入れると発散が収まったことから、渦粘性が効いて流速が馴染んだのかもしれません。
ということは、流体の物性値の粘性係数を大きくしても同じことが起こりますし、メッシュサイズを小さくすることで発散が改善できるかもしれません。
また、浮力による自然対流で乱流ができているであれば浮力を発生させる重力を０にすると計算が進むとも考えられます。
ということでこれらを考慮して設定を変えてスタディしました。
まず結論からまとめると以下のようになりました。
〇：設定する、×：設定しない
結果は予想通りです。
チュートリアルをコピー
今回はブシネスク近似を使わないソルバ「buoyantSimpleFoam」を使います。空気などで温度差が15℃以上ある場合は、ブシネスク近似での誤差が1%以上出てくるので温度差がある解析にはブシネスク近似は使わないようにしましょう。
以下のコマンドで熱流体解析ができるチュートリアルをコピーします。
cp -r $FOAM_TUTORIALS/heatTransfer/buoyantSimpleFoam/hotRadiationRoom base環境変数として「$FOAM_TUTORIALS = /opt/OpenFOAM/OpenFOAM-v2012/tutorials」が割り当てられています。
フォルダ構成はこんな感じ。
.
├── 0
│ ├── G
│ ├── T
│ ├── U
│ ├── alphat
│ ├── epsilon
│ ├── k
│ ├── nut
│ ├── p
│ └── p_rgh
├── Allclean
├── Allrun
├── constant
│ ├── boundaryRadiationProperties
│ ├── g
│ ├── radiationProperties
│ ├── thermophysicalProperties
│ └── turbulenceProperties
└── system
├── blockMeshDict
├── controlDict
├── fvSchemes
└── fvSolutionパターンを変えるためにケースファイルをコピーします
cp -r base case1
cp -r base case2
cp -r base case3
cp -r base case4
cp -r base case5
cp -r base case6
cp -r base case7輻射の設定
輻射の設定はオフにする
「constant/radiationProperties」は輻射を設定するときに使うファイルなので、今回は以下のようにオフにしておきます。
constant/radiationProperties
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: v2012 |
| \\ / A nd | Website: www.openfoam.com |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
location "constant";
object radiationProperties;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
radiation on;
radiationModel none;
/*
radiationModel P1;
// Number of flow iterations per radiation iteration
solverFreq 1;
absorptionEmissionModel constantAbsorptionEmission;
constantAbsorptionEmissionCoeffs
{
absorptivity absorptivity [0 -1 0 0 0 0 0] 0.5;
emissivity emissivity [0 -1 0 0 0 0 0] 0.5;
E E [1 -1 -3 0 0 0 0] 0;
}
scatterModel none;
sootModel none;
/*
// ************************************************************************* //
変更した個所は以下です。
radiationModel none;乱流モデルの設定
乱流モデルは「constant/turbulenceProperties」で行います。
RANSのk-εモデルの場合
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: v2012 |
| \\ / A nd | Website: www.openfoam.com |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
location "constant";
object turbulenceProperties;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
simulationType RAS;
RAS
{
RASModel kEpsilon;
turbulence on;
printCoeffs on;
}
// ************************************************************************* //
乱流モデル無しの場合
simulationType laminar;
/*
simulationType RAS;
RAS
{
RASModel kEpsilon;
turbulence on;
printCoeffs on;
}
*/
乱流モデルを使わない場合は、デフォルトの設定をコメントアウトして、simulationType をlaminarとします。
重力の設定
重力の設定は「constant/g」で行います。
constant/g
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: v2012 |
| \\ / A nd | Website: www.openfoam.com |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class uniformDimensionedVectorField;
location "constant";
object g;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [0 1 -2 0 0 0 0];
value (0 0 9.8);
// ************************************************************************* //
デフォルト乗せていだと重力はzの負の方向に与えられています。
重力無しの場合はvalueの値を変更します。
dimensions [0 1 -2 0 0 0 0];
value (0 0 0);粘性係数の設定
流体の物性値の設定は「constant/thermophysicalProperties」で行います。
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: v2012 |
| \\ / A nd | Website: www.openfoam.com |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
location "constant";
object thermophysicalProperties;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
thermoType
{
type heRhoThermo;
mixture pureMixture;
transport const;
thermo hConst;
equationOfState perfectGas;
specie specie;
energy sensibleEnthalpy;
}
pRef 100000;
mixture
{
specie
{
molWeight 28.9;
}
thermodynamics
{
Cp 1000;
Hf 0;
}
transport
{
mu 1.8e-05;
Pr 0.7;
}
}
// ************************************************************************* //
粘性係数大きくする場合
粘性係数を1000倍にしてみます。
transport
{
mu 1.8e-2;
Pr 0.7;
}メッシュサイズを1/2倍
メッシュ生成は「system/blockMeshDict」で行います。
こちらはヘキサメッシュでメッシュを作るOpenFOAMのユーティリティのblockMeshです。
メッシュサイズを変更したい場合は以下の部分で行います。
blocks
(
hex ( 0 1 5 4 16 17 21 20) ( 5 5 5) simpleGrading (1 1 1)
hex ( 1 2 6 5 17 18 22 21) (10 5 5) simpleGrading (1 1 1)
hex ( 2 3 7 6 18 19 23 22) (80 5 5) simpleGrading (1 1 1)
hex ( 4 5 9 8 20 21 25 24) ( 5 10 5) simpleGrading (1 1 1)
hex ( 6 7 11 10 22 23 27 26) (80 10 5) simpleGrading (1 1 1)
hex ( 8 9 13 12 24 25 29 28) ( 5 40 5) simpleGrading (1 1 1)
hex ( 9 10 14 13 25 26 30 29) (10 40 5) simpleGrading (1 1 1)
hex (10 11 15 14 26 27 31 30) (80 40 5) simpleGrading (1 1 1)
hex (16 17 21 20 32 33 37 36) ( 5 5 15) simpleGrading (1 1 1)
hex (17 18 22 21 33 34 38 37) (10 5 15) simpleGrading (1 1 1)
hex (18 19 23 22 34 35 39 38) (80 5 15) simpleGrading (1 1 1)
hex (20 21 25 24 36 37 41 40) ( 5 10 15) simpleGrading (1 1 1)
hex (21 22 26 25 37 38 42 41) (10 10 15) simpleGrading (1 1 1)
hex (22 23 27 26 38 39 43 42) (80 10 15) simpleGrading (1 1 1)
hex (24 25 29 28 40 41 45 44) ( 5 40 15) simpleGrading (1 1 1)
hex (25 26 30 29 41 42 46 45) (10 40 15) simpleGrading (1 1 1)
hex (26 27 31 30 42 43 47 46) (80 40 15) simpleGrading (1 1 1)
);例えば「 hex ( 0 1 5 4 16 17 21 20) ( 5 5 5) simpleGrading (1 1 1)」の( 5 5 5)を(10 10 10)と倍すれば分割数が倍になります。
分割数を倍にしたい場合は、以下のようにします。
blocks
(
hex ( 0 1 5 4 16 17 21 20) (10 10 10) simpleGrading (1 1 1)
hex ( 1 2 6 5 17 18 22 21) (20 10 10) simpleGrading (1 1 1)
hex ( 2 3 7 6 18 19 23 22) (160 10 10) simpleGrading (1 1 1)
hex ( 4 5 9 8 20 21 25 24) (10 20 10) simpleGrading (1 1 1)
hex ( 6 7 11 10 22 23 27 26) (160 20 10) simpleGrading (1 1 1)
hex ( 8 9 13 12 24 25 29 28) (10 80 10) simpleGrading (1 1 1)
hex ( 9 10 14 13 25 26 30 29) (20 80 10) simpleGrading (1 1 1)
hex (10 11 15 14 26 27 31 30) (160 80 10) simpleGrading (1 1 1)
hex (16 17 21 20 32 33 37 36) (10 10 30) simpleGrading (1 1 1)
hex (17 18 22 21 33 34 38 37) (20 10 30) simpleGrading (1 1 1)
hex (18 19 23 22 34 35 39 38) (160 10 30) simpleGrading (1 1 1)
hex (20 21 25 24 36 37 41 40) (10 20 30) simpleGrading (1 1 1)
hex (21 22 26 25 37 38 42 41) (20 20 30) simpleGrading (1 1 1)
hex (22 23 27 26 38 39 43 42) (160 20 30) simpleGrading (1 1 1)
hex (24 25 29 28 40 41 45 44) (10 80 30) simpleGrading (1 1 1)
hex (25 26 30 29 41 42 46 45) (20 80 30) simpleGrading (1 1 1)
hex (26 27 31 30 42 43 47 46) (160 80 30) simpleGrading (1 1 1)
);これで再度blockMeshを実行すると分割数がx,y,z方向で倍になっています。
メッシュ数は8倍(2×2×2)です。
計算の実行
以上のように設定をひとつひとつ見てきましたが、設定を変えて実行してみましょう。
計算がうまくいったCase4で進めます。
blockMesh
buoyantSimpleFoamメッシュはこのようになっています。
断面を切って適当に温度スケールを変えるとこんな感じ。
実験データとの比較
いずれにしても適当に乱流モデルを設定して上手く行きましたが、いつまでも付きまとうのが「解析の妥当性あるの？」という疑問です。
あくまでCAE解析は設計の支援ツールであるため、実現象を再現してこそ意味を持ちます。
実験データをとることが難しい場合はネットで公開されている実験データと解析を比較するのも良いでしょう。
解析したいモデルが実物の形状が良いという場合は、自身で実験をする必要がありますが、今回のように熱流体は「空調の流れは？床は断熱なのか？」など外部要因が多いので、実験と全く同じ条件で解析するのは不可能かもしれません。
しかし、解析でモデル化できないにしても空間の温度分布の測定などの実験条件を定めておくことは実機と解析に違いが出た場合のヒントとしては非常に重要な情報となります。
OpenFOAMをはじめるなら
OpenFOAMをはじめるなら手にしておきたい日本語の書籍がこちらです。
OpenFOAMをインストールしてこれから深く学んでいこうという方や、使い慣れてもレベルアップのために何回も読みたくなるくらい参考になる内容です。
初心者の方にはチュートリアルを動かしながら学んでいくのが良いでしょう。こちらはOpenFOAMをインストールしてチュートリアルを少し触っていくところからはじめて最後に重合格子をやるなど面白い内容です。
Twitter➡@t_kun_kamakiri
ブログ➡宇宙に入ったカマキリ（物理ブログ）
youtube➡理系ブログ
Instagram➡kamakiri1225