タグ: blockMesh

SwiftツールによるOpenFOAM®用メッシュ作成

個人事業主となって、公開でギャラの貰える仕事としては初仕事になる表題の講習会が終了しました。

オープンCAE学会でこれまでに実施してきた講習会の実績からすると、予想以上の集客になってしまい、しかも内容が内容だけに、補助できるサポート要員もない講習会で、一部の人には消化不良になってしまったかもしれません。しかし、Twitterなどでの反応を見る限り、大方の人には及第点をつけて頂けたのでないかと、まずはほっとしております。

一番の反省点としては、Macユーザーさん対応でした。ライブDVDでやると決めた時点から、Macユーザーさんの事はスッパリ抜けていましたが、直前になって「Macユーザーはどうしたらいいんだ？」という問い合わせがあり、仮想マシンでなんとかなるでしょう・・・と、その先のことは全く考えてなかったもんで。ここはいっそ、「Macは駄目です」と回答しておけば、講習会としてはよほどスイスイやれたでしょう・・・ということにはなります。

まぁ、今回、Macは仮想マシンでやってもらって、一番の問題「Macのマウスにはホイールが無い」点に関する配慮さえやっておけば、さほど大きな問題もなくやれそうだ・・・ということが判ったということで、ここは前向きに考えております。

で、そこんところは、どうすんだ？　・・・　を帰ってから、色々調べたんですが、blenderをMacで使う方法に関する情報があまりに少ないことには、少々驚きました。それでもまぁ、色々調べたところ、多分ですが、ホイールを回す相当の機能は、テンキーを使ってほぼ実現できるので、Macの場合もその方法（テンキー）で出来るんじゃないかと思いますが、何せMac使いではないので、Macのテンキーがどうなっているかも知らないので、後はMacユーザーさんにおまかせするしかありません。

一応、マウス操作に関する説明文中に、ホイール操作と同等のテンキー操作方法について内容を追加したものを、講習テキストの電子版（ｐｄｆ）として、ここに公開しておきます。⇒SwiftツールによるOpenFOAM®用メッシュ作成

SwiftBlock作成例：四角柱周りの流れ

公開コンサル案件にて、四角柱周りの流れを計算したいとの依頼があったので、このところ注力中のSwiftBlockを使ってメッシュ作成に挑戦してみました。

caseファイルと使用法の説明

  1. caseファイルは上のように展開されるはず（２Dの場合）
  2. pillar2DRSwift.blend をblenderで開く（3Dの場合はpillar3DSwift.block）
  3. SwiftBlockのメニューからwriteボタンを押してblockMeshDictを作成する
  4. blockMeshDictの出力先は、constant/polyMesh の下
  5. （メッシュの分割数やグレーディングを変更したい場合は、上記blockMeshDictを直接編集して下さい）
  6. blockMeshの実行
  7. （0フォルダの下のFiledファイルはSwiftBlockにて設定したpatch名と整合済。境界条件は便宜的に設定しただけなので、流入条件などを変更したければ、直接編集にて変更して下さい）
  8. （constant および system フォルダの下にある各種設定ファイルは、tutorials/incompressible/pimpleFoam/Tjunction にあったものをそのまま流用しています。適宜変更して使って下さい）
  9. pimpleFoamの実行
  10. SwiftBlockによるpatch定義をする前の状態のブロックモデルも同梱してあるので、SwiftBlockの設定を変更したい場合には、こちらから始めて下さい。

ブロック作成例(2D)

  • サイズは適当に作っただけですが、この程度の作業であれば10分もあれば十分に可能。

メッシュ作成例（2D)

  • メッシュ数

nPoints: 3160

nCells: 1468

計算結果例（2D）

ブロック作成例(3D)

メッシュ作成例（3D)

  • メッシュ数

nPoints: 48040

nCells: 43632

計算結果例（3D）

補足

本稿を見て、ブロックの作成方法に「何でこの形？」という疑問を持たれた方も居ると思います。たとえば、２Dの場合だと、以下のような形

のほうが、よほど簡単でチュートリアルにふさわしいと思われるしょう。これは実際にこのモデル（pillar2D.blend）も同梱してあるので試してもらうのが良いと思うし、現実に最初に試したのはこのモデルでした。しかしどうしてもメッシュがうまく作れませんでした。うまく作れないの意味は、中央の四角柱部分までメッシュが出来てしまうという事です。

まぁ、後はblockMeshDictを手修正でやれないこともなかったんですが、これはSwiftBlockのバグかもしれないと思い、作者へフィードバックすべきかどうか考えるうちに、これはバグではなく、仕様であるとの結論に至りました。

つまり、SwiftBlockがblockMeshDictを生成する原理は、線画の構造体の中から六面体構造を検出して、その部分をblockMeshDictの作法に則ってblock定義出力しているにすぎない、という事です。（あくまで推測ですが・・・）

なので、四角形柱部分を六面体構造で定義してあるので、block定義されてメッシュが作られてしまうという事です。

そこで、四角柱部分をデータ構造としては六面体構造でない形で定義する、２Dの場合五角形柱、3Dの場合には六角形柱になるようにしてブロック図を作成してみたら、予想通りうまくいった・・・という訳です。

3次元のブロック構造を作るに際しても、面の引き伸ばしでブロックを拡張していくと、内部格子点間での線情報が欠落する場合があり、そうなるとメッシュが出来ないということになります。メッシュが予想通りに切れない場合は、本観点でもって作成したブロックの線情報を見直す、というのが最優先の着眼点でしょう。

caseファイル

SwiftBlock

表題のツールを使うと、blender上で図形を見ながらGUIでblockMeshDictを作成できるという嬉しさがあって、速報ではごく簡単な直方体だけの場合の作成例を紹介した。

しかし、wikiのページに掲載された例では、円柱用のblockMeshDictも作成できるとあり、更にはそこからダウンロードできるexampleケースを使うと、下図のような、そこそこに複雑な六面体のブロック集合モデル（この種のモデルであれば、blenderで結構簡単に作成できる）から、

作成したblockMeshDictを使ってメッシュ作成すると・・・

何と、見事な六面体メッシュが出来てくれたので、新鮮な驚きでありました。

２つのオブジェクトを使用してblockMeshDictを作成

そこで、どうやったらこれが出来るかのかを調べたところ、まず、このblenderのモデルには、２つのオブジェクトが存在しているという点。

モデルを開いた時の状態は、blockモデルしか表示されていないが、画面右上のオブジェクトツリーを見ると、もう一つpartモデルがあることがわかり、このeyeマークをクリックするとモデル画面上に表示される。

要するに、blockモデルを使ってblockMeshDictを作成するんですが、その際に曲線部などはpartモデルにフィッティングするよう、blockMeshDict中のedge定義をpoyLineで自動作成してくれているんだということです。

２つのオブジェクトで節点共有させておくのがミソ

そこまでのところは、すぐに判ったのですが、問題は、それでどうやってモデルを作ったら良いのか？　という点です。マニュアルによれば、Set edgesにチェックマークが入っており、

If an edge in your block structure has both vertices co-located with vertices from the typed object, SwiftBlock will find out the closest path in the mesh of the object

blockモデルと、partモデルでco-located（節点共有？）させなさい・・・ということらしい。で、co-locateさせるにはどうしたら良いのかというと、

To co-locate a block vertex with a vertex of a Geometry object (to allow for non-straight block edges), use Blender’s snapping magnet in vertex mode.

という説明でおしまい。これも、blender使いには、こういう説明で判るかもしれないが、凡人には無理。あれこれ調べて、ようやく判明しました。

 節点共有させる方法

円柱用のblockMeshを作成する場合を例として説明する。

まず、円柱と直方体を別々のオブジェクトとして作成する。上下端面の高さは一致させてあり、直方体は、円柱の内部に収まる形で配置。

円柱用のblockMeshは、中央直方体の４つの側面に台形ブロックを継ぎ足す方法が普通のやり方。最終的に下図のような、５つのブロックがあれば、円柱用のblockMeshDictを作成できることになる。

これを実現するやり方として、普通に考えつくのは、側面を伸長させる方法だろう。Edit Modeの面編集モード(Face select mode)にて、

ひとつの側面を選択状態にして、Eキーを押す（又は、Meshメニューから、Extrude Region を選択する）と、マウスの動きと共に面の法線方向へブロックが伸長する。

適当な長さまで引き伸ばす。マウス左ボタンをクリックすれば確定する。ついで、頂点編集モードにて、

適当な長さにまで拡張したブロックの端点を選択した状態の上図にて、問題はこの点を円柱モデル上の点、たとえば下図

の選択点と、どうやって一致させるかという事。ここでポイントになるのが、中央の選択点の右にある赤白の破線の丸印で示されるカーソルです。これを使って選択点を所望の点と一致させます。

Meshメニューから、Snapを選択すると・・・

selectionをどうしたい、cursorをどうしたい、という選択肢が出てきます。ここまでくればおわかりでしょう・・・

まずは、カーソルを選択点（円柱のメッシュ点）に持っていきます（下図）。

この時点で、円柱を対象のEdit Mode になっていますが、Tabキーを押して、Object Mode にして、対象Objectを円柱からブロックに変更し、もう一度Tabキーを押して、再度Edit Modeに戻す。

そうすると、ブロックのEdit状態であるので、選択点はカーソルから離れて、ブロック上のどこかへ移動する。ここで、共有したい節点を選択し、再度、MeshメニューからSnapを選択

今度は、Selection to Cursor を選択してやれば、


上図のようになって、ブロックの角端点が、円柱のメッシュ点と一致してくれるようになる訳です。残りの３点も同様のやり方で円柱のメッシュ点と一致させる。さらに中央直方体の残り３つの側面も同じように拡張して、最終的に上のほうで見た５つのブロックが出来れば完成ということです。

 

回転円筒容器内の水面挙動（その１）

OpenFOAMユーザーズグループで投稿のあった、回転円筒容器内の水面挙動に関し、MRFInterFOAMを使って計算する方法について調べてみた。⇒確かに、なかなか一筋縄で計算できるものでなかったようで、本稿（その１）では、まずは計算してみた篇です。

普通に予想される結果（放物面状の液面）とは異なっていますが、一概にまるでおかしい・・・とも言い切れない面があると思います。どうでしょう？少し偏心して回転すると、これに似た挙動をしてくれそうな気もします（実は、ビール造りで、大きなバケツ状タンクに投入した酵母イーストを攪拌する際によく経験しています）。

最後のほうで、何故こうなったのか、どうすれば予想結果が得られそうかの考察をしていますが、まずはここまで、普通にやればこうなったというところを、以下具体的にどうやって計算したかの解説です。

なお、あちらではsnappyHexMeshでメッシュ作成していましたが、こちらではblockMeshのみでやりました。（当初、提示されていたモデルではメッシュが辛そうだったので）

チュートリアルケースの変更概要

基本的なやり方は、チュートリアルケース（tutorials/multiphase/MRFInterFoam/mixerVessel2D）をそのままコピーしてケース名を変更、一部のファイルを変更するという作業です。ファイル別に、実施した修正作業の全体概要を以下に示しておきます。

以下、修正内容の要点を記しておきます。

メッシュ作成

メッシュ作成には、拙作のブロックメッシュ作成マクロを使いました。
マクロファイル（*.m4)をblockMeshDictのあるpolyMeshフォルダと同じ場所に置いておけば、ケースのルートフォルダにて、以下のコマンドを実行すれば、メッシュの出来上がりです。

$ m4 < constant/polyMesh/blockMeshDict_Pipe-4.m4 > constant/polyMesh/blockMeshDict
$ blockMesh

マクロのパラメタ設定部分

メッシュの出来上がりイメージ

境界定義

あちらでは、上面を開口にしていましたが、ここでは（面倒なので）全周壁面としての計算にした。多分、あんまり結果は変わらないんではないかと・・・

上の方法で作成したメッシュの境界（constant/polyMesh/boundary）は、全面がdefaultFaces（type empty）となっているので、

これを、walls (type wall) に変更する。手修正も簡単だが、メッシュを変更の都度やり直すのは面倒だし、自動化もやりたい。ここはcreatePatchコマンドでやることにした。

systemフォルダの下に、以下の内容のcreatePatchDictを仕込んでおけば、

以下のコマンド実行で、一発変換してくれます。

$ createPatch

境界条件

流速

圧力

液体比率

回転領域の設定

MRF（Multi Reference Frame）系のソルバーを使うには、回転領域に相当する部分のセルにVolume名をつけておく必要がある。色々なやり方があるが、ここではtopoSetコマンドを用いる。

systemフォルダの下に、以下の内容のtopoSetDictを仕込んでおけば、

$ topoSet

$ setsToZone -noFlipMap

というコマンドにて、全領域に、rotate_area というVolume名が付けられることになります。

 

初期水面高さの設定

水面の初期状態はsetFieldコマンドにて設定します。systemフォルダの下に、以下のsetFieldDict を仕込んでおけば、

$ cp 0/alpha1.org 0/alpha1

$ setField

にて、ｚ方向高さ１メートルの部分までが、水（alpha1=1）で満たされたことになる。

回転条件の設定

回転条件は、constantフォルダの下の、MRFZonesファイルで記述しておく

原点（0 0 0）を通る、ｚ軸（0 0 1）回りに、5.25 rad/s ということです。

 

計算制御（controlDict）

endTimeが400と非常に大きい値にしてあります。チュートリアルケースでは４でした。なので、途中でリスタート計算できるよう、

sttartFrom   latestTime;

と変更しました。

writeIntervalは、0.1くらいにしておかないと、アニメーション表示に耐えられません（0.1secで30degほどの回転角になってしまうので）。400secの計算ともなると、4000コマの出力になってしまいます。ここはご自分の環境にあわせて書き換えて下さい。

deltaTは1e-3となっていますが、あくまで初期値です。以下、

にて、自動時間刻みで計算は進行します。この部分の値は、チュートリアルケースの値を変更していません。

最後に、function probe を追加しています。

底面の上方0.1メートルの中心部（0.0 0.0 0.1）と外周部（0.45 0.0 0.1）における流速(U)と圧力(p_rgh)を全計算ステップにてモニター出力します。

これがないと、計算が準定常状態になったかどうかの判断が出来なかったからです。

計算実行

添付のAllrunを実行するだけです。

#!/bin/sh

m4 < constant/polyMesh/blockMeshDict_Pipe-4.m4 > constant/polyMesh/blockMeshDict
blockMesh
createPatch -overwrite
topoSet
setsToZones -noFlipMap
cp 0/alpha1.org 0/alpha1
setFields
MRFInterFoam | tee log.MRFInterFoam

# —————————————————————– end-of-file

計算結果の概要

上段に前述のprobe点における圧力履歴をプロットし（緑色が外周、赤色が中央）、下段に任意の時刻における液面を表示しています。

このように、回転が始まってから30sec（約25回転）まで、ほとんど変化らしいものは見られず、35secくらいで、ようやくさざ波のようなものが現れて、それが次第に大きくなる。300secを過ぎたあたりでようやく準安定状態といってよさそうです。

probe点圧力履歴を、350-400sec間だけプロットし直したのが下の図です。

この図から読み取って、最後の２周期分、約8sec間0.1sec毎のデータで作成したのが、冒頭に掲載したアニメーションです。

ちなみに計算時間は単体計算で約6時間

供試マシンは、ちょっと前の記事にて紹介したLinuxマシン

CPU:Core-i7 950
Memory:24GB
OS:Ubuntu-10.04

上に構築したMint12ベースの仮想マシン（ubuntu-11.10 相当）でやりました。

何はともあれ、お試しあれ

右（⇒）フレームのDownLoadカラムの　rotatingCylinder　です。

考察

当初、通常のファンなどで実施するMRF計算の経験から、数回転から１０回転分の計算をすれば出来るだろう・・・と思って、endTimeをチュートリアルと同じ（4sec）か、数倍（20sec）まで見ておけばよかろうと始めたのですが、まるで液面が変化しませんでした。メッシュや計算パラメタに誤りがあるのではないかと、あれこれ変更確認し、随分回り道をしましたが、結局計算時間が足りなかったということでした。

通常やっている内部にファンを配したような計算では、内部物体表面の移動境界条件によって、流体の回転方向対流が生じることになりますが、本ケースの場合は、壁面の摩擦だけで流体を回転させることになるので、流体全体の回転運動が生じるようになるまでは相当の時間が必要になるということのようです。

ただ、ようやく計算できたとはいえ、冒頭に記したように予想した液面が放物面状にはなってくれませんでした。この結果が必ずしも不合理ではないかもしれないとしても、何故こういうことになったのかを考えると、本現象のレイノルズ数的なもの（慣性力/粘性力の比を代表する指標）を見ていく必要がありました。

たとえば、代表寸法を、円筒の直径（D=1）、代表速度を周速（U=RΩ,R=0.5,Ω=5.25）として考えると、

nu nu [ 0 2 -1 0 0 0 0 ] 1e-06;

なので、

Re=UD/nu
=2R^2Ω/nu
=0.5 x 10e6

と、非常に大きな慣性力支配の現象になっているということがわかります。なので粘性力による回転流れが出る前に慣性力による振動が出てしまう。そういう流れを層流計算している・・・という点にも無理があったということかもしれません。

したがって、予想した放物面の液面形状を得るには、以下２つのアプローチが考えられます。（他にもあるかもしれませんが・・・）

  1. 低レイノルズ数計算
  2. 乱流計算

以下、（その２）へ続く。（近日公開予定）