snappyHexMeshを用いた角柱周りの流体解析について

126 views
Skip to first unread message

S.TAKEHARA

unread,
Apr 12, 2024, 3:41:02 AMApr 12
to OpenFOAM
お世話になっております。
大学の研究でOpenFOAMを用いての流体解析を行おうとしています、初学者のTAKEHARAと申します。

まずは円柱周りの流れについてカルマン渦を確認しようと考え、
以下のサイトを参考に、非定常解析ソルバーとしてpimpleFoamを用いて円柱周りの流れ解析を行い、正常に結果を得ることができました。

そこで、次に比較のために、角柱周りの流れ解析を行おうと考えています。
上記サイトとでは、blockMeshのみでメッシュを生成しているようなのですが、
それは難解だと感じたため、snappyHexMeshを用いてメッシュを作成しています。
snappyHexMeshを用いることで角柱メッシュを生成後、pimpleFoamを行うのですが、計算が発散してしまっている状況です。

以下メッシュ情報です。
形状0.06 m×0.13 m(縦×横)
一辺2.0×10^-3 mの角柱を左端から0.03m位置に配置
z軸方向の壁面についてはemptyとし、
x軸方向の壁面はpatch
y軸方向の壁面と角柱表面はwallとしています。

メッシュ生成後、checkMeshも行っていますが、そこでは
Mesh OK.
となっております。

初学者のため、境界条件や乱流モデルの設定等のミスがあるかもしれません。
お忙しい中恐縮ですが、エラー原因のほど、ご教授のほどいただけると幸いです。

pimpleFoamによるエラーメッセージ:
[stack trace]
=============
#1  Foam::sigFpe::sigHandler(int) in /Volumes/OpenFOAM/OpenFOAM-v2306/platforms/darwin64ClangDPInt32Opt/lib/libOpenFOAM.dylib
#2  _sigtramp in /usr/lib/system/libsystem_platform.dylib
#3  Foam::DICPreconditioner::calcReciprocalD(Foam::Field<double>&, Foam::lduMatrix const&) in /Volumes/OpenFOAM/OpenFOAM-v2306/platforms/darwin64ClangDPInt32Opt/lib/libOpenFOAM.dylib
#4  Foam::DICPreconditioner::DICPreconditioner(Foam::lduMatrix::solver const&, Foam::dictionary const&) in /Volumes/OpenFOAM/OpenFOAM-v2306/platforms/darwin64ClangDPInt32Opt/lib/libOpenFOAM.dylib
#5  Foam::lduMatrix::preconditioner::addsymMatrixConstructorToTable<Foam::DICPreconditioner>::New(Foam::lduMatrix::solver const&, Foam::dictionary const&) in /Volumes/OpenFOAM/OpenFOAM-v2306/platforms/darwin64ClangDPInt32Opt/lib/libOpenFOAM.dylib
#6  Foam::lduMatrix::preconditioner::New(Foam::lduMatrix::solver const&, Foam::dictionary const&) in /Volumes/OpenFOAM/OpenFOAM-v2306/platforms/darwin64ClangDPInt32Opt/lib/libOpenFOAM.dylib
#7  Foam::PCG::scalarSolve(Foam::Field<double>&, Foam::Field<double> const&, unsigned char) const in /Volumes/OpenFOAM/OpenFOAM-v2306/platforms/darwin64ClangDPInt32Opt/lib/libOpenFOAM.dylib
#8  Foam::PCG::solve(Foam::Field<double>&, Foam::Field<double> const&, unsigned char) const in /Volumes/OpenFOAM/OpenFOAM-v2306/platforms/darwin64ClangDPInt32Opt/lib/libOpenFOAM.dylib
#9  Foam::fvMatrix<double>::solveSegregated(Foam::dictionary const&) in /Volumes/OpenFOAM/OpenFOAM-v2306/platforms/darwin64ClangDPInt32Opt/lib/libfiniteVolume.dylib
#10  Foam::fvMatrix<double>::solveSegregatedOrCoupled(Foam::dictionary const&) in /Volumes/OpenFOAM/OpenFOAM-v2306/platforms/darwin64ClangDPInt32Opt/lib/libfiniteVolume.dylib
#11  main in /Volumes/OpenFOAM/OpenFOAM-v2306/platforms/darwin64ClangDPInt32Opt/bin/pimpleFoam
#12  start in /usr/lib/dyld
=============
zsh: floating point exception  pimpleFoam
prism.zip

Hideaki Kominami

unread,
Apr 12, 2024, 7:53:50 AMApr 12
to open...@googlegroups.com
TAKEHARAさん

kominamiです。

自分の環境はOpenFOAMのv2206で、TAKEHARAさんのOF10とは違いますが、とりあえず見てみました。

./Allrunスクリプト内にsnappyHexMeshコマンドがありません。TAKEHARAさんがコマンドを手打ちしていると解釈して、自分も手打ちしました。
自分の環境のv2206では、snappyHexMeshのaddLayersControlsのところで、minMedialAxisAngleの設定が無いというエラーメッセージがでるのですが・・・

>メッシュ生成後、checkMeshも行っていますが、そこでは
>Mesh OK.
>となっております。

snappyHexMeshのところでは、エラーが出ていませんか?
snappyHexMeshの動作でエラーがあっても、その後のcheckMeshが通ることがあります。

以上、よろしくお願いします。

2024年4月12日(金) 16:41 S.TAKEHARA <souta...@gmail.com>:
--
このメールは Google グループのグループ「OpenFOAM」に登録しているユーザーに送られています。
このグループから退会し、グループからのメールの配信を停止するには openfoam+u...@googlegroups.com にメールを送信してください。
このディスカッションをウェブ上で閲覧するには https://groups.google.com/d/msgid/openfoam/e98ceb59-5394-48de-8639-4dd007493118n%40googlegroups.com にアクセスしてください。

mmer547

unread,
Apr 12, 2024, 9:46:18 AMApr 12
to OpenFOAM
komioamiさんの補足ですが、snappHexMeshを並列で実行している場合、やり方によっては
ベースのblockMeshに対してcheckMeshを実行している場合があるので、
その辺教えてもらえるとエラー原因検討できるかと思います。
2024年4月12日金曜日 20:53:50 UTC+9 kominami:

S.TAKEHARA

unread,
Apr 13, 2024, 2:43:25 AMApr 13
to OpenFOAM
kominamiさん

お忙しい中、動作確認ありがとうございます。
環境についての記載が漏れていました。OpenFOAM v2306でコマンドは手打ちしております。
blockMesh→snappyHexMesh→decomposePar→mpirun -np 8 pimpleFoam -parallel
で動かしております。

>snappyHexMeshのところでは、エラーが出ていませんか?

当初、minMedialAxisAngleを90で設定していたのですが、 minMedialAxisAngleの設定が無いというエラーメッセージが表示されました。
また、エラーが表示されていてもparaviewで確認するとメッシュ自体の生成はできているように見えました。
ですので、このエラーについては解消法が分からないことも相まって、放置してしまっている状況です。

2024年4月12日金曜日 22:46:18 UTC+9 mmer547:

S.TAKEHARA

unread,
Apr 13, 2024, 2:46:48 AMApr 13
to OpenFOAM
mmer547さん

snappHexMeshは並列せずに実行しています。
pimpleFoamについては並列で実行予定です。

2024年4月13日土曜日 15:43:25 UTC+9 S.TAKEHARA:

nakagawa

unread,
Apr 13, 2024, 5:00:07 AMApr 13
to OpenFOAM
TAKEHARAさん

なかがわです。

まずは,作成されたメッシュを目視で確認することをお勧めします。

角柱まわりのメッシュが添付図のようになっています。角部付近のひし形に近いセルは,発散の要因になりやすいです。着目するであろう物体の対称性とメッシュの対称性の違いがあるなど,あとから結果の解釈などに苦労するようにも感じます。

cylinder_mesh_closeup.png


基本的な形状であれば,blockMeshでメッシュを使う方がよい場合が多いです。snappyHexMeshDictをいじり倒してメッシュを調整していると,blockMeshで計画的に作った方が早かったとなることがあります。いろいろ勉強していくと,急がば回れという言葉の重みを感じることがあります。OpenFOAMを使い始めたときは,どちらも経験することで,メッシュ作成力が向上するのですが。

blockMeshで正方形柱を対象とするとき,複数ブロックを使う方法(参考にされて円柱と同様)があります。それ以外に,単一ブロックでメッシュを作成した後に,不要部分(角柱)を抜き取る(topoSetとsubsetMeshを使う)方法があります。一見,面倒に感じる方法ですが,使えるようになると便利です。

オープンCAE勉強会@富山の第82回の講習会テキストが参考になると思います。下記からご覧ください。

角柱まわりのメッシュを細分化したり,レイヤーを追加したりすることが,OpenFOAMのユーティリティでできました。
2013年8月31日 第13回 オープンCAE勉強会@富山 の講習資料と例題ケースを確認してはどうでしょうか。上記リンクからアクセスできます。

計算が発散するときは,結果を保存する時間を短く設定して,発散直前の結果を可視化することをお勧めします。圧力や速度の大きな部分を可視化したり,thresholdなどで抽出して可視化し,メッシュとの関係を把握すると対策が立てやすいです。

がんばってください。

なかがわ




2024年4月13日土曜日 15:46:48 UTC+9 S.TAKEHARA:

S.TAKEHARA

unread,
Apr 13, 2024, 5:37:39 AMApr 13
to OpenFOAM
なかがわさん

お忙しい中、ご丁寧に返信くださりありがとうございます。

>角柱まわりのメッシュが添付図のようになっています。角部付近のひし形に近いセルは,発散の要因になりやすいです。着目するであろう物体の対称性とメッシュの対称性の違いがあるなど,あとから結果の解釈などに苦労するようにも感じます。

解析条件ももちろん関係はあると思いますが、それよりもやはりメッシュの形状が発散の原因となっている可能性が高いでしょうか。

>基本的な形状であれば,blockMeshでメッシュを使う方がよい場合が多いです。snappyHexMeshDictをいじり倒してメッシュを調整していると,blockMeshで計画的に作った方が早かったとなることがあります。いろいろ勉強していくと,急がば回れという言葉の重みを感じることがあります。OpenFOAMを使い始めたときは,どちらも経験することで,メッシュ作成力が向上するのですが。

>blockMeshで正方形柱を対象とするとき,複数ブロックを使う方法(参考にされて円柱と同様)があります。それ以外に,単一ブロックでメッシュを作成した後に,不要部分(角柱)を抜き取る(topoSetとsubsetMeshを使う)方法があります。一見,面倒に感じる方法ですが,使えるようになると便利です。

そうなのですね。今回のケースのような角柱だとかなり基本的な形状になると思うので、ご教授いただいた複数ブロックを使う方法と不要部分を抜き取る方法について詳しく調べてみようと思います。

>2013年8月31日 第13回 オープンCAE勉強会@富山 の講習資料と例題ケースを確認してはどうでしょうか。上記リンクからアクセスできます。

講義資料については以前より何度もお世話になっていますが、再度確認して例題ケースから自分の環境で試してみるようにしていきます。

今回スレッドに初投稿させていただいたのですが、今後も利用させていただこうと考えていますので、何卒よろしくお願いいたします。
2024年4月13日土曜日 18:00:07 UTC+9 nakagawa:

k k

unread,
Apr 14, 2024, 1:57:03 AMApr 14
to open...@googlegroups.com
解決済みでしたらすみません。
snappyHexMeshでメッシュ作成して、自分の手元でもやってみました。
※ただし、v2312を使っていますのでOF10とは異なりますが、設定はできるだけ同じにしました。

【やったこと】
・2次元解析と思ったのでblockMeshのzの大きさを以下に変更。
計算には意味ないと思いますが、z方向に長すぎて見にくいため。
vertices
(
    (-0.03 -0.03 -0.005)
    ( 0.1 -0.03 -0.005)
    ( 0.1  0.03 -0.005)
    (-0.03  0.03 -0.005)
    (-0.03 -0.03 0.005)
    ( 0.1 -0.03 0.005)
    ( 0.1  0.03 0.005)
    (-0.03  0.03 0.005)
);
同時にsnappyHexMeshDictの以下を変更
 locationInMesh ( 0.05  0.02 0); // Inside point
・メッシュが細かすぎて時間がかかるので、以下に変更
blocks
(
    hex (0 1 2 3 4 5 6 7) (100 100 1) simpleGrading (1 1 1)
);

・snappyHexMeshの再分割と境界層の相性が良くなく、境界層が上手く入りにくいので以下に変更
    featureAngle 130;//30;//60;
ただし、OF10は比較的境界層がきれいに入ったと思うので、あまり気にしなくても良いかもしれません。

これで計算エラーはなくなります。
=============
上記までの設定は、blockMeshでメッシュを粗めに切ったため、渦が再現できていない(ani_3.gif)のでメッシュを細かくする必要があります。
計算時間を考えると、すべてを領域を細かくする必要はないです。
角柱背後の渦を再現するためには領域分割で対応すれば文献値と近くなると思います。
例えば、領域のメッシュ再分割はsnappyHexMeshDictのこちらを設定します。
    refinementRegions
    {
    }


snappyHexMeshを使った角柱の2次元解析であれば、こちらの方がやりたいことに近いかと思います。

根本のエラーの原因は、わかりませんが、
・snappyHexMeshでメッシュ生成できていないのに、pimpleFoamを実行していたから
・メッシュ品質が悪い
くらいしか思いつきません。

以上はsnappyHexMeshの例ですが、blockMeshのみでもできると思うので色々試してみてください。

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

2024年4月13日(土) 18:37 S.TAKEHARA <souta...@gmail.com>:
このディスカッションをウェブ上で閲覧するには https://groups.google.com/d/msgid/openfoam/f905ed83-8a7a-46cf-913a-3e6bc447981an%40googlegroups.com にアクセスしてください。
prism.png
メッシュ.png
ani_3.gif

Hideaki Kominami

unread,
Apr 14, 2024, 6:30:33 AMApr 14
to open...@googlegroups.com
TAKEHARAさん

kominamiです。

【1】

>環境についての記載が漏れていました。OpenFOAM v2306でコマンドは手打ちしております。

>blockMesh→snappyHexMesh→decomposePar→mpirun -np 8 pimpleFoam -parallel
>で動かしております。

snappyHexMesh -overwrite というオプションが必要です。
-overwriteオプションが無いと、controlDictのtimeStepで指定された時間刻みでメッシュが作成されます。snappyHexMeshの動作チェックのときは各段階の途中経過を残すために-overwriteオプション無しで使う場合もありますが、ソルバー計算を行うときには-overwriteオプションが必須です。


【2】
controlDictのstartFromが時刻0からスタートする設定になっている場合は、TAKEHARAさんの入力コマンドだと、blockMesh直後のメッシュから計算されます。したがって計算エラーの発生にはsnappyHexMeshは関係ありません。(今回は、こちらの設定ですね。)

controlDictのstartFromが最終時刻からスタートする設定になっている場合は、なかがわさんが添付した図のように、発散しやすいメッシュを使ってることになります。

【3】
kkさんが言及していますが、2次元モデル(奥行方向のメッシュ分割数が1)なのに、奥行方向の寸法が非常に長いのが気になります。アスペクト比の大きなメッシュは計算の安定性を損ないます。今回は2次元モデルだから、checkMeshでエラーが出ず、計算も問題ないという可能性はあるでしょう。(自分は、2次元モデルで奥行の寸法をやたら大きくしないためで実際のところは分かりません。)

以上、よろしくお願いいたします。


2024年4月14日(日) 14:57 k k <kusuko...@gmail.com>:
このディスカッションをウェブ上で閲覧するには https://groups.google.com/d/msgid/openfoam/CAEAawWEwFSksB8QgzK11MsTjt9%2BUG7nDkU-bHBy%3DALfksbph4Q%40mail.gmail.com にアクセスしてください。

S.TAKEHARA

unread,
Apr 15, 2024, 1:44:34 AMApr 15
to OpenFOAM
kkさん

エラーの解消まで確認していただき、ありがとうございます。
ご教授いただいた設定を元に私の環境でも実行してみたところ、無事に計算をすることができました。
(checkMeshではメッシュのチェックに失敗しているものの、計算自体は実行できているようです。checkMeshについては以下に添付しておきます。)

>上記までの設定は、blockMeshでメッシュを粗めに切ったため、渦が再現できていない(ani_3.gif)のでメッシュを細かくする必要があります。

渦についてはおそらくレイノルズ数が達していなかったと思われます。試しに流速を0.1 m/s →1.0 m/sに変更してみたところ渦が確認できました。
メッシュについてはご教授いただいた通り、refinementRegionsを変更してみようと思います。

スクリーンショット 2024-04-15 13.39.47.png

お陰様で計算ができるようになり、本当に助かりました。
ご指導いただき、改めて御礼申し上げます



checkMesh(hex (0 1 2 3 4 5 6 7) (200 200 1) simpleGrading (1 1 1) で設定しています)
--------------------------------------------------------------
Create time

Create mesh for time = constant

Check mesh...

Time = constant

Mesh stats
    points:           83988
    faces:            167960
    internal faces:   85744
    cells:            42224
    faces per cell:   6.00853
    boundary patches: 6
    point zones:      1
    face zones:       0
    cell zones:       0

Overall number of cells of each type:
    hexahedra:     42094
    prisms:        0
    wedges:        0
    pyramids:      0
    tet wedges:    0
    tetrahedra:    0
    polyhedra:     130
    Breakdown of polyhedra by number of faces:
        faces   number of cells
            6   15
            9   110
           12   5

Checking topology...
    Boundary definition OK.
 ***Total number of faces on empty patches is not divisible by the number of cells in the mesh. Hence this mesh is not 1D or 2D.
    Cell to face addressing OK.
    Point usage OK.
    Upper triangular ordering OK.
    Face vertices OK.
    Number of regions: 1 (OK).

Checking patch topology for multiply connected surfaces...
                   Patch    Faces   Points  Surface topology
                   y_min      200      402  ok (non-closed singly connected)
                   y_max      200      402  ok (non-closed singly connected)
                   x_min      200      402  ok (non-closed singly connected)
                   x_max      200      402  ok (non-closed singly connected)
            defaultFaces    81104    82056  ok (non-closed singly connected)
                   prism      312      390  ok (non-closed singly connected)
                    ".*"    82216    82290      ok (closed singly connected)


Checking faceZone topology for multiply connected surfaces...
    No faceZones found.

Checking basic cellZone addressing...
    No cellZones found.

Checking basic pointZone addressing...
               PointZone  PointsBoundingBox
            frozenPoints       0(1e+150 1e+150 1e+150) (-1e+150 -1e+150 -1e+150)

Checking geometry...
    Overall domain bounding box (-0.03 -0.03 -0.005) (0.1 0.03 0.005)
    Mesh has 2 geometric (non-empty/wedge) directions (1 1 0)
    Mesh has 2 solution (non-empty) directions (1 1 0)
 ***Number of edges not aligned with or perpendicular to non-empty directions: 515
  <<Writing 871 points on non-aligned edges to set nonAlignedEdges
    Boundary openness (-1.46618e-17 -7.33791e-17 -7.77721e-14) OK.
    Max cell openness = 3.31625e-16 OK.
    Max aspect ratio = 15.1515 OK.
    Minimum face area = 7.13537e-10. Maximum face area = 6.88249e-06.  Face area magnitudes OK.
    Min volume = 1.78385e-12. Max volume = 2.09807e-09.  Total volume = 7.79601e-05.  Cell volumes OK.
    Mesh non-orthogonality Max: 85.2877 average: 5.69704
   *Number of severely non-orthogonal (> 70 degrees) faces: 480.
    Non-orthogonality check OK.
  <<Writing 480 non-orthogonal faces to set nonOrthoFaces
    Face pyramids OK.
    Max skewness = 1.68308 OK.
    Coupled point location match (average 0) OK.

Failed 1 mesh checks.

End
2024年4月14日日曜日 19:30:33 UTC+9 Hideaki Kominami:

S.TAKEHARA

unread,
Apr 15, 2024, 2:05:41 AMApr 15
to OpenFOAM
kominamiさん

今までは最新のメッシュをコピペして利用していたのですが、-overwriteオプションをつけることでそれが不要になるのですね。
早速取り入れてみます。

奥行方向の寸法については円柱解析を行った際に参考にさせていただいたサイトを元に設定していました。
初学者のため、あまり勝手に変更するべきではないと思っていたのですが、ご教授いただいたように、確かに二次元モデルの場合、長くすることに意味はないですね。
アスペクト比が大きいメッシュだと安定性が損なわれてしまうというのも新しい学びでとても勉強になりました。
kkさんにも指摘いただいたように、z方向を短くすることで計算が流れたので、今回はその要因もあったのではないかと推測しております。

kominamiさんにもとても多くの新しい学びをいただき、感謝しております。
本当にありがとうございます。
2024年4月15日月曜日 14:44:34 UTC+9 S.TAKEHARA:
Reply all
Reply to author
Forward
0 new messages