【OpenFOAM TIPS】objectRegistryにオブジェクトを登録する
はじめに
よく忘れるので備忘録としてまとめます。
OpenFOAMの計算実行時には、objectRegistryオブジェクトがメモリ内のオブジェクトを制御しています。ここに登録してあるものには色んな場所からアクセスることができます。この記事では、coded境界条件の初回呼び出し時にフィールドを登録し、以降は同一のそれを呼び出せるようにする、という事をやっていきます。
計算環境
- OS : Ubuntu 20.04 (WSL / Windows 10 Home)
- OpenFOAM : v2012
メリット
以下のようなメリットが考えられます。
- 一つのオブジェクトを色々な場所から共有できる
- 出力のタイミングを他のオブジェクトと合わせることができる
- 場として出力できればParaViewで他のフィールドと同様に可視化できる
このやり方はcoded等の書き捨てコードに限らず、例えば独自の境界条件を作る場合などにも有効です。例えば速度と圧力の境界条件間で何かオブジェクトを渡したい場合に、計算を制御するobjectRegistryにそのオブジェクトを登録することでメモリ内で実現できたりします。
例として用いたケース
ここでは例として、inletの質量流量をvolScalarFieldクラス"massInlet"の境界値として各時刻で保存する場合を考えます。
以下のような脈動する水と油の二相流を例に説明します。（※本記事ではこの計算自体には触れません。）
なおここで使用したケースはGithub(github.com/inabower/objectRegistryTest)で公開しています。
脈動を表す境界条件
Uでは以下のようにcodedFixedValueで脈動を設定しているとします。内容としては時間によってx方向の速度が脈動するような条件です。上下の位置で分布を持たせたりもしています。
boundaryField
{
inlet
{
type codedFixedValue;
name wave;
value uniform (0 0 0);
code
#{
const fvMesh& mesh = patch().boundaryMesh().mesh();
const Time& t = mesh.time();
vectorField ans(patch().size(), Zero);
forAll(ans, fi)
{
ans[fi].x() = 5.0* (1.0-cos(t.value()*3.141592*2.0))*(0.25-pow(patch().Cf()[fi].y(), 2.0));
}
operator==(ans);
#};
}
...(他の境界)
}
ここのcodedの中にいろいろ追加していきます。具体的にはmassInletというフィールドを作成し、その境界値に密度と速度を掛けた値を追加することを考えます。
ダメな例
通常のcreateFields.Hのようにこの中でvolScalarFieldを定義することをやってみます。
boundaryField
{
inlet
{
type codedFixedValue;
name wave;
value uniform (0 0 0);
code
#{
const fvMesh& mesh = patch().boundaryMesh().mesh();
const Time& t = mesh.time();
vectorField ans(patch().size(), Zero);
forAll(ans, fi)
{
ans[fi].x() = 5.0* (1.0-cos(t.value()*3.141592*2.0))*(0.25-pow(patch().Cf()[fi].y(), 2.0));
}
operator==(ans);
const word massInletName = "massInlet";
volScalarField massInlet
(
IOobject
(
massInletName,
mesh.time().timeName(),
mesh,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
mesh,
dimensionedScalar(dimless, Zero)
);
scalarField& massInletRef = massInlet.boundaryFieldRef()[patch().index()];
const scalarField& alpha = patch().lookupPatchField<volScalarField, scalar>("alpha.water");;
forAll(massInletRef, fi)
{
massInletRef[fi] = (alpha[fi] * 1000 + (1.0 - alpha[fi]) * 800) * ans[fi].x();
}
#};
}
...(他の境界)
}
これだと一応計算は回りますが、massInletは結果としては出力されません。これはvolScalarField オブジェクトの寿命がこの関数内で終了していることに由来します。オブジェクトを保存して後のイタレーションでも使用するためにはobjectRegistryに登録する必要があります。
よさげな例
以下のように、ポインタをobjectRegistryに保存すると、以降の計算でも同一のオブジェクトを呼び出すことができ、さらにcontrolDictなどで指定したwriteFileのルールに従い他のフィールドと同じタイミングで出力してくれます。
boundaryField
{
inlet
{
type codedFixedValue;
name wave;
value uniform (0 0 0);
code
#{
const fvMesh& mesh = patch().boundaryMesh().mesh();
const Time& t = mesh.time();
vectorField ans(patch().size(), Zero);
forAll(ans, fi)
{
ans[fi].x() = 5.0* (1.0-cos(t.value()*3.141592*2.0))*(0.25-pow(patch().Cf()[fi].y(), 2.0));
}
operator==(ans);
fvMesh& meshRef = const_cast<fvMesh&>(mesh);
const word massInletName = "massInlet";
if(!meshRef.foundObject<volScalarField>(massInletName))
{
autoPtr<volScalarField> massInletPtr(
new volScalarField
(
IOobject
(
massInletName,
meshRef.time().timeName(),
meshRef,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
meshRef,
dimensionedScalar(dimless, Zero)
)
);
meshRef.objectRegistry::store(massInletPtr);
}
const scalarField& massInletFld = patch().lookupPatchField<volScalarField, scalar>(massInletName);
scalarField& massInletRef = const_cast<scalarField&>(massInletFld);
const scalarField& alpha = patch().lookupPatchField<volScalarField, scalar>("alpha.water");;
forAll(massInletRef, fi)
{
massInletRef[fi] = (alpha[fi] * 1000 + (1.0 - alpha[fi]) * 800) * ans[fi].x();
}
#};
}
...(他の境界)
}
ここでは以下の動作を行っています。
- objectRegistryに"massInlet"が登録されていない場合
- volScalarFieldクラスの"massInlet"を新しく生成しそのポインタを
massInletPtrとして取り出す massInletPtrポインタを引数にしてobjectRegistryに登録する
- volScalarFieldクラスの"massInlet"を新しく生成しそのポインタを
patch().lookupPatchField関数を用いて、objectRegistryに登録されているmassInletのpatchの値を呼び出す- 内容を変更する
注意すべき点をいくつか書いておきます。
① 発動するタイミング
if(!meshRef.foundObject<volScalarField>(massInletName))
objectRegistryに登録するのは必ず一回にしなくてはいけません。そのためここではfound関数を用いて、すでに登録されていない場合のみstoreするようにしています。
② store関数
autoPtr<volScalarField> massInletPtr( new volScalarField(io) );
meshRef.objectRegistry::store(massInletPtr);
登録するためのobjectRegistryクラスのstore関数を用います。ここではfvMeshがobjectRegistryクラスを継承しており、OpenFOAMの計算の根幹のobjectRegistryを呼び出すことができるため、meshRef.objectRegistry::store(massInletPtr);のようにしています。
ちなみにautoPtr<>はOpenFOAMで用意しているスマートポインタです。可能な限りC++ネイティブのものではなくこれを使うようにした方が個人的には思っています。（OpenFOAMの作者の意図に反したエラーを回避できるため）
③ const_cast
const fvMesh& mesh = patch().boundaryMesh().mesh(); fvMesh& meshRef = const_cast<fvMesh&>(mesh);
ここ(fvPatchField)からは以下のようにしてfvMeshクラスのポインタにアクセスできます。
- const fvPatch& patch()
- const fvBoundaryMesh& patch().boundaryMesh()
- const fvMesh& patch().boundaryMesh().mesh()
ただこれらはconstであるため、const_castで内容を変更可能にします。
そのほか応用例
空のオブジェクトを登録する操作を入れることで、最初の一回かどうかのフラグとして使用することもできます。 ここでは最初の１回目のイタレーションのみ速度0、それ以外は速度1にしています。
boundaryField
{
inlet
{
type codedFixedValue;
name wave;
value uniform (0 0 0);
code
#{
vector ans(Zero);
const fvMesh& mesh = patch().boundaryMesh().mesh();
fvMesh& meshRef = const_cast<fvMesh&>(mesh);
bool isFirstTime = false; //ここ
if(!mesh.foundObject<IOdictionary>("firstIterFlag_"+patch().name()))
{
autoPtr<IOdictionary> flagPtr(
new IOdictionary
(
IOobject
(
"firstIterFlag_"+patch().name(),
meshRef.time().constant(),
meshRef,
IOobject::NO_READ,
IOobject::NO_WRITE
)
)
);
meshRef.objectRegistry::store(flagPtr);
isFirstTime = true; //ここ
}
// 最初のイタレーションでない場合はx方向の速度を1.0
if(!isFirstTime) ans.x() = 1.0;
operator==(ans);
#};
}
...(他の境界)
}
計算の安定化などで役立つこともあるかもしれません。
最後に
もしも何かの役に立ったら幸いです。
参考
- 公式Doxygen
snappyとヘキサメッシュ
はじめに
OpenFOAMではsnappyHexMeshというユーティリティがあり、これを用いてメッシングを行う人も多いかと思います。
そしてその多くはblockMeshで作成した構造格子に対して、計算したいモデルのstlファイルを用いてsnappyHexMeshを行っているかと思います。
ただこのsnappyHexMeshはヘキサメッシュであれば何にでも適応可能であり、事前にSALOMEなどでヘキサメッシュを用意できるのであれば結構大きなメリットがあります。
今回はこの例を紹介したいと思います。
作業環境
OS : Ubuntu 18.04 LTS OpenFOAM : v1912 SALOME : 9.4.0
アルミ缶の外部流れ
あまりないかもしれませんが、アルミ缶の外部流れを計算したい場合を考えます。
この場合の一般的な戦略としては、以下のようになるかと思います。
- アルミ缶の形状のSTLファイルを作る
- 外部流れ場を
blockMeshで作成 snappyHexMeshでアルミ缶部分を切り抜く
このとき、境界層がうまく作れなかったりアルミ缶周りのメッシュがあまり綺麗でなかったり局所的に小さなメッシュができたりなどの苦労をすることが想定されます。
そして苦労する割には肝心のアルミ缶周りのメッシュはポリヘドラやプリズムが混ざったものになってしまいます。そこで今回は以下のような手順を考えます。
- SALOMEでアルミ缶周りを重視した円筒形のヘキサメッシュを作成する
- 床面のSTLファイルを作成する
snappyHexMeshで床より下を切り抜いたメッシュを作成する
こうすることで最も重要なアルミ缶周りの綺麗なメッシュを作成することができました。 目的に合わせて上面や横面も同時にカットすることもできます。
変な形の撹拌翼
あまりないかもしれませんが以下のような変な形の撹拌翼のメッシュを作りたい場合を考えます。
この場合もblockMesh+形状STLよりも、SALOMEなどで外側ヘキサメッシュを作成した後に複雑な部分のみをSTLにしてsnappyHexMeshすると以下のように綺麗に作成することができます。
最後に
是非お試しください。
SALOMEでOpenFOAM用のwedgeメッシュを作る
はじめに
OpenFOAMでは軸対象のモデルを計算するwedgeというモードが備わっています。 このメッシュの作成にちょっとコツがいるためまとめてみました。
作業環境
- OS : Ubuntu18.04 LTS
- SALOME : 9.4.0
- OpenFOAM : v1912
作成したメッシュ
以下のような350mLアルミ缶のモデルについて行いました。
途中までは以下の記事で作成したものを使用します。
このSHAPERモジュールで作成してGEOMETRYモジュールにExportしたShell_1を元に作成していきます。
wedge
wedgeモデルは軸対象のモデルに対して軸回転方向を無視したものです。 軸を上から見て10°以内の扇型のメッシュについて計算を行うことができます。 この角度が大きすぎるとセル中心の位置が面中心の位置とずれてしまい、小さすぎるとセル体積が小さくなり誤差の原因になります。 私はなんとなく5°(1/72)で使用することが多いです。
GEOMETRY
1. 中心部分をPartition
SHAPERなどで作成した四角形分割済みの面について、 Point→Lineで中心部分1メッシュ分の線を引き、この線でShellをPartitionします。
2. 扇形にする
Revolution：5°で扇形にし、その後にRotation：-2.5°で位置を中心に合わせます。
3. 中心部分だけグループ分け
必然的にプリズム型になる中心部分は他のヘキサメッシュとは異なるアルゴリズムを使用するため名前をつけておきます。 Extrusionというベース面のメッシュをそのまま押し出すようなアルゴリズムを使用するため、Extrusionする立体とそのベースとなる面をグループ化しておきます。
4. 線のグループ分け
同じ分割数の線をそれぞれグループ化します。 Operations → Blocks → Propagateを使うと便利です。
MESH
5. 線だけメッシュ
まずは線だけメッシングします。 そのために先程の扇形オブジェクトに対してCreate Meshを行い、3Dと2DのアルゴリズムはNoneにしておきます。 そしてCreate SubMeshで4で作成したグループ全てに対して分割数を指定します。 なお扇形の円周方向と中心部の外内方向のメッシュは分割数1（分割なし）にします。
6. base面を単一三角形としてメッシング
分割数1の辺×3の三角形であるbase面を単一の三角形にします。
”メッシングしない”というオプションがあれば楽なのですが（見つけてないだけかも）、ここではNETGEN2Dアルゴリズムを使用します。
ここでAllow Quadrangleにチェックを入れるとなぜか単一三角形として生成してくれます。
7. 単一三角形を押し出し
3Dアルゴリズムの中の3D Extrusionを先程のextに適用します。
この状態でComputeを行うことで単一三角形を押し出したプリズムメッシュでext領域をメッシングできます。
8. その他の領域をヘキサメッシュ
他の部分についてはヘキサメッシュを行います。
最初にアルゴリズムを指定せずに作った大元のメッシュに対してEdit meshを行います。
2DアルゴリズムをQuadrangle Mapping、3DアルゴリズムをHexahidron(i,j,k)にすることで残りの部分をヘキサメッシュにできます。
9. GEOMETRYに戻って面のグループ分け
OpenFOAMでpatchになる面に名前をつけていきます。
私の場合はいつもwedge_Aとwedge_Bにしていますが、これらはそれぞれ別のグループになっているのであれば何でもOKです。
10. MESHモジュールに戻って読み込み→出力
メッシュを選択して右クリックメニューのCreate Groups from Geometryで先程のグループをメッシュに読み込ませます。
その後、右クリックメニューのExport→UNV fileでUNV形式にエクスポートします。
今回はwedge.unvというファイルを出力しました。
11. OpenFOAM形式に変換
OpenFOAMの適当なチュートリアルをコピーしてきてその中に先程のwedge.unvをコピーし、以下の操作を行います。
ideasUnvToFoam wedge.unv
renumberMesh -overwrite
ideasUnvToFoamはUNV形式のメッシュをOpenFOAM形式に変換するユーティリティです。
renumberMeshはメッシュの順番を整えるユーティリティで、SALOMEで出力されるメッシュはメッシュの順番がバラバラな場合が多いためこれを行うことで計算速度が早くなったりします。
また、constant/polyMesh/boundaryの中身を以下のように編集します。
/*--------------------------------*- C++ -*----------------------------------*\ | ========= | | | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox | | \\ / O peration | Version: v1912 | | \\ / A nd | Website: www.openfoam.com | | \\/ M anipulation | | \*---------------------------------------------------------------------------*/ FoamFile { version 2.0; format ascii; class polyBoundaryMesh; location "constant/polyMesh"; object boundary; } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // 6 ( wedge_A { type wedge; nFaces 1571; startFace 3055; } wedge_B { type wedge; nFaces 1571; startFace 4626; } ...
これでwedge面として指定できました。
12. メッシュの確認
ParaViewとcheckMeshで確認したメッシュは以下のようになっていました。
Checking geometry...
Overall domain bounding box (-1.113691e-18 -0.0014394398 -2.0157582e-12) (0.032968591 0.0014394398 0.117)
Mesh has 2 geometric (non-empty/wedge) directions (1 0 1)
Mesh has 3 solution (non-empty) directions (1 1 1)
Wedge wedge_A with angle 2.5 degrees
Wedge wedge_B with angle 2.5 degrees
All edges aligned with or perpendicular to non-empty directions.
Boundary openness (1.1379343e-16 1.9791924e-16 5.5112448e-19) OK.
Max cell openness = 2.3610339e-16 OK.
Max aspect ratio = 4.4778269 OK.
Minimum face area = 4.3577871e-08. Maximum face area = 6.6817001e-06. Face area magnitudes OK.
Min volume = 2.9051911e-11. Max volume = 7.2329335e-09. Total volume = 5.2055678e-06. Cell volumes OK.
Mesh non-orthogonality Max: 56.766778 average: 23.621304
Non-orthogonality check OK.
Face pyramids OK.
Max skewness = 1.3492111 OK.
Coupled point location match (average 0) OK.
いい感じですね。
最後に
wedgeにするとメッシュ数が格段に減るため軸回転方向以外の解像度を大幅に高めることができます。 ただし軸回転方向の影響を無視できない場合には結果が大きく変わるためご注意ください。
SALOMEのSHAPERモジュールでヘキサメッシングしてみた
はじめに
先日SALOMEの9.4.0がリリースされました。 前回のSALOME 9.3.0から追加されたSHAPERモジュールが、今回の9.4.0でバグ修正や機能追加が行われており結構使いやすくなってきています。 今回はこのSHAPERを用いてヘキサメッシュを行う方法について紹介します。
作業環境
- OS : Ubuntu 18.04 LTS
- SALOME : 9.4.0
作成したメッシュ
こんなアルミ缶のメッシュを作成しました。
- メッシュ数：104768
- 形状：ヘキサメッシュ
- 最大非直交性：56.77
- 平均非直交性：23.14
- 最大歪： 1.349
- 最大アスペクト比：5.921
SHAPERモジュール
SHAPERモジュールとは
元々SALOMEでは作図のGEOMETRYモジュールとメッシングのMESHモジュールが搭載されていました。 GEOMETRYモジュールでは点→線→面→立体と拡張していきながら3D図形を作成していきますが、作成前に座標などを正確に把握する必要がありました。 一方で今回のSHAPERモジュールはスケッチを起点にインタラクティブな作図を行うことができるようになります。 将来的にはGEOMETRYモジュールがSHAPERに取り込まれるとのウワサもあります。
公式ドキュメントはこちら
公式チュートリアル動画もあります。
Salome Shaper Tutorial 1 : Modeling a Parametric Nut
ボタンの概要
このボタンがSHAPERです。
操作の感覚としてはFreeCADに近いかと思います。 スケッチ面を決定してその中の図形を、長さや角度などで拘束しながら作図していきます。
具体的には実例を見せながら説明していきます。
1. スケッチ面を作成して画面を合わせる
まずスケッチ面を作成します。 今回はX-Z面を使用します。 ただこの指定により作成されるスケッチ面は自動的に1m×1mの大きさになり、今回作成するアルミ缶の大きさに対して大きすぎるので、直線を一本引いてその長さを指定し、それに合わせて拡大していき画面のスケールを合わせます。
2. 線を引く
次に大体の線を引いていきます。
点や線の上に重なった状態で点を打つとそこで拘束されます。
もし外れてしまった場合は拘束しなおせばOKです。
また水平や垂直はある程度自動で設定されます。
もし水平垂直が不要な場合にはVやHを選んでDELキーで消去します。
3. 線を拘束する
長さや角度などを拘束していくことにより図形を目的のものに近づけていきます。 完全に拘束されて自由度が0になると線が緑色になり、この状態を目標にしていきます。 拘束が過剰になるとそれ以上変更できなくなるので、元に戻すか該当の拘束を削除します。
4. 分割用の内部線を作成する
今回はヘキサメッシュを目的としていますので、面は四角形で分割する必要があります。 この四角形分割でSHAPERモジュールが力を発揮します。 面をきれいな四角形に分割するために、例えば水平拘束をしたり、少しだけ角度をつけたり、同じ長さにしたりなど、調整をしながら線を描いていくことができます。 もともとのGEOMETRYモジュールにおける作図では予め座標の計算などを行わなくてはいけなかったところをインタラクティブに行うことができるようになりました。
5. Partsを作成してGEOMETRYモジュールに送る
スケッチからオブジェクトに起こすにはPartsという機能を用います。
まずPart → New Partで新しいPartスペースを作成して、そこに今回はスケッチからシェルオブジェクトを作成します。
作成した後はExport to GEOMでGeometryモジュールに送ります。
GEOMETRY
ここからは元々の機能ですのでサクサクやっていきます。
こちらの詳細は過去の記事で紹介しています。
6. 分割用の内部面を作成する
回転した円筒を六面体で分割できるように8角柱ベースで分割しています。
7. 内部線のグループわけ
Operations → Blocks → Propagateで同じ分割数にすべき線を自動でグループ化してくれます。
少なくともこのグループで線の分割数を指定してあげればヘキサメッシュを切ることができます。
今回は更にこのグループ同士Unionすることで後のMESH工程での指定を簡単にしています。
MESH
8. 線毎に分割数を指定してメッシュ作成
普通にやります。
できました
最後に
なんとなくSHAPERで綺麗なメッシュが切れるのが分かりましたでしょうか。 まだ出たての新機能であまり資料はありませんが是非試してみてください。
OpenFOAM-v1906のインストール〜ParaView5.7.0とPython3を添えて〜
はじめに
先日OpenFOAM-v1906がリリースされました。せっかくなので最新のParaViewと一緒にインストールする方法をまとめてみました。
注意
- この記事の内容（コマンド等）に関して一切の責任を負いかねます。
sudoを用いたコマンドが頻出します。入力ミス、コピペミス、私のミスなどで重大なエラーが生じる可能性がありますのでコマンド実行の際にはその意味を確かめてから実行してください。- インストール時の環境変数がAnacondaと競合することがあります。念の為OpenFOAMインストール時にはAnacondaは無効化するのが無難です。
- 【未解決】本記事では以下の機能をインストールできませんでした。
- ParaViewの
-mpiオプション：7/12現在ParaViewで並列計算を有効にするとインストールできない事例を確認しています。（参考） - PVFoamReader：ParaViewのOpenFOAMケース表示プラグインであるPVFoamReaderをインストールできませんでした。
- ParaViewの
環境
- OS : Ubuntu 18.04 LTS（VirtualBox内）
- OpenFOAM : v1906
- ParaView : 5.7.0-RC1
- Python : 3.6.8
- インストール日 : 2019/07/11〜7/13
参考
以下のサイトを見ながらインストールしました。困ったら全部良く読みましょう。
- OpenFOAM® Installation From Source
- doc/Build.md · master · Development / OpenFOAM-plus · GitLab
- Installation/Linux/OpenFOAM-v1806/Ubuntu - OpenFOAMWiki
- はじめに
- 1. Requirements
- 2. OpenFOAMなどのダウンロード
- 3. ParaView-5.7.0-RC1をダウンロード
- 4. ParaViewのバージョンを指定するファイル
- 5. 環境変数の設定
- 6. ParaViewのコンパイル
- 7. OpenFOAMのコンパイル
- 8. インストールの確認
- 最後に
1. Requirements
まずはOpenFOAMやParaViewを動かすためのライブラリなどをインストールします。(Python3も入れます)
sudo apt update && sudo apt -y upgrade sudo apt install build-essential flex bison cmake zlib1g-dev libboost-system-dev libboost-thread-dev \ libopenmpi-dev openmpi-bin gnuplot libreadline-dev libncurses-dev libxt-dev \ qt5-default libqt5x11extras5-dev libqt5help5 qtdeclarative5-dev qttools5-dev \ libqtwebkit-dev freeglut3-dev libqt5opengl5-dev texinfo \ libscotch-dev libcgal-dev python python-dev \ libglu1-mesa-dev python3 python3-dev libglu1-mesa-dev qt4-dev-tools
仮想マシンの場合
VirtualBoxなどの仮想マシンの場合は以下の操作が必要になるそうです。（OpenFOAM wiki）
sudo apt install mesa-utils scons llvm-dev sudo sed -i~orig -e 's/# deb-src/deb-src/' /etc/apt/sources.list sudo apt update sudo apt build-dep mesa
2. OpenFOAMなどのダウンロード
各ソースをダウンロード＆解凍していきます。
# ディレクトリの作成 mkdir OpenFOAM && cd OpenFOAM # OpenFOAM-v1906のダウンロードと解凍 wget https://sourceforge.net/projects/openfoam/files/v1906/OpenFOAM-v1906.tgz tar -xzf OpenFOAM-v1906.tgz # ThirdParty-v1906のダウンロードと解凍 wget https://sourceforge.net/projects/openfoam/files/v1906/ThirdParty-v1906.tgz tar -xzf ThirdParty-v1906.tgz # decomposeParのためのmetisがはいっていないためダウンロードと解凍 cd ThirdParty-v1906 wget http://glaros.dtc.umn.edu/gkhome/fetch/sw/metis/metis-5.1.0.tar.gz tar -xzf metis-5.1.0.tar.gz
3. ParaView-5.7.0-RC1をダウンロード
公式サイトから"ParaView-v5.7.0-RC1.tar.gz"ダウンロードします。
そしてこれをさきほど作成された~/OpenFOAM/ThirdParty-v1906の中で解凍し、解凍されたディレクトリの名前をParaView-v5.7.0に変更します。
そうするとディレクトリの構成はmetisも含めて以下のようになっているかと思います。
4. ParaViewのバージョンを指定するファイル
現在ThirdParty-v1906ではデフォルトでParaView-v5.6.0をインストールする設定になっています。これを変更するにはOpenFOAM側の設定を変更する必要があります。 例えばターミナルに以下のコマンドを入力しファイルを開きます。
gedit ~/OpenFOAM/OpenFOAM-v1906/etc/config.sh/paraview
ここで57行目のParaView_VERSION=を5.7.0に変更することでParaViewのバージョンを指定することができます。
... # # Note # When _foamAddLib is unset (eg, called from makeParaView or from foamPV): # - the ParaView_VERSION variable is retained. #------------------------------------------------------------------------------ # USER EDITABLE PART: Changes made here may be lost with the next upgrade ParaView_VERSION=5.7.0 ParaView_QT=qt-system cmake_version=cmake-system ...
※後に出てくる./makeParaViewの際に./makeParaView 5.7.0のように引数にバージョンを記述することで同様に5.7.0をインストールすることもできます。
しかしこれだとOpenFOAM側のutilitiesであるparaFoamで呼び出すことができません。
今回はこれも踏まえてOpenFOAM側のconfigファイルを変更しました。
5. 環境変数の設定
ターミナル(bash)起動時にOpenFOAM1906の環境変数を読み込むためにbashrcに追記を行います。
echo "alias of1906='source \$HOME/OpenFOAM/OpenFOAM-v1906/etc/bashrc'" >> $HOME/.bashrc echo "of1906" >> $HOME/.bashrc . $HOME/.bashrc
No completion added for /home/inabower/OpenFOAM/OpenFOAM-v1906/platforms/linux64GccDPInt32Opt/bin
... incorrect platform, or not yet compiled?
上のようなメッセージが出たら準備完了です。更に確認したい場合にはgedit $HOME/.bashrcなどのコマンドで.bashrcファイルを開いて最後の行に
alias of1906='source $HOME/OpenFOAM/OpenFOAM-v1906/etc/bashrc' of1906
が追加されていればOKです。aliasによりof1906と入力すればOpenFOAM-v1906の環境をセットできるようになります。そのためv1906をデフォルトで使用しないのであればof1906での行は削除して構いません。
その場合はターミナルを開くたびにof1906入力することになります。
ちょっといじる
私の環境だとこのまま./makeParaViewすると2時間後くらいに以下のエラーが出てしまいコンパイルができませんでした。
CMake Error at /home/inabower/OpenFOAM/ThirdParty-v1906/ParaView-v5.7.0/CMake/ParaViewClient.cmake:522 (message): Failed to convert servermanager XML: xmlpatterns: could not exec '/usr/lib/qt5/bin/xmlpatterns': No such file or directory make[2]: *** [Plugins/AcceleratedAlgorithms/paraview_help/AcceleratedAlgorithms_doc.xslt] Error 1 make[1]: *** [Plugins/AcceleratedAlgorithms/CMakeFiles/AcceleratedAlgorithms_doc.dir/all] Error 2 make: *** [all] Error 2 Command exited with non-zero status 2
苦肉の策ですが、以下のようにqt4から持ってくるとエラーは解消しました。
sudo cp /usr/lib/x86_64-linux-gnu/qt4/bin/xmlpatterns /usr/lib/qt5/bin/
※ もし/usr/lib/qt5/bin/xmlpatternsが既にあるようでしたらこの操作は問題ありません。
- こことかを見ながらqt5のxmlpatternsをaptからインストールしようと試みたのですが叶いませんでした。もし方法をご存じの方がいましたら教えてください。
6. ParaViewのコンパイル
新しいターミナルを開き、以下のコマンドでParaViewをインストールしていきます。-python -python-lib -python-includeではPython3のライブラリの場所を指定します。
「ファイルが見つからない」と言われた場合にはfind /usr/lib -name "libpython3*"などのコマンドで探してみましょう。
cd $WM_THIRD_PARTY_DIR ./makeParaView -python -python-lib /usr/lib/x86_64-linux-gnu/libpython3.6m.so.1.0 -python-include /usr/include/python3.6m/ > log.makeParaView
終了まで数時間かかります。なお新しいターミナルで以下のコマンドを入力することでインストールの進行状況を見ることができます。
tail -f $WM_THIRD_PARTY_DIR/log.makeParaView
...
-- Installing: /home/inabower/OpenFOAM/ThirdParty-v1906/platforms/linux64Gcc/ParaView-5.7.0/lib/python3.6/site-packages/paraview/vtk.py
====
Installation complete for paraview-5.7.0 with qt-5.9.5
ParaView_DIR=/home/inabower/OpenFOAM/ThirdParty-v1906/platforms/linux64Gcc/ParaView-5.7.0
You may need to update the OpenFOAM environment by running:
wmRefresh
====
Done
というようなメッセージができたらインストール成功です。
7. OpenFOAMのコンパイル
いよいよOpenFOAMをコンパイルしていきます。
foam ./Allwmake -j 4
4はコンパイルの並列数です。各自のCPUの数に合わせて数字を入れてください。
2019-07-13 12:32:07 +0900 ======================================== Finished compile of cfMesh with OpenFOAM-v1906 Gcc system compiler linux64GccDPInt32Opt, with SYSTEMOPENMPI openmpi-system 2019-07-13 12:32:07 +0900 ======================================== OpenFOAM-v1906 Gcc system compiler linux64GccDPInt32Opt, with SYSTEMOPENMPI openmpi-system api = 1906 patch = 0 bin = 285 entries lib = 133 entries ========================================
こんなのが表示されたら完了です。（途中でインストールを中断したりしたので数などが異なる可能性があります。）
8. インストールの確認
ちゃんとインストールされているかどうかを確認していきます。
軽いソルバー(potentialFoam) + paraFoam
今回インストールしたものの中で、まず軽いソルバーの実行とその結果の表示ができるかを確認します。
ソルバーはpotentialFoamで試してみて、それをparaFoamで表示します。
paraFoamはParaViewがインストールされた状態でOpenFOAMをインストールすると搭載される機能で、OpenFOAMのケースをコマンド一つで開くことができるようになります。
mkdir $WM_PROJECT_USER_DIR && cd $WM_PROJECT_USER_DIR cp -r $FOAM_TUTORIALS/basic/potentialFoam/pitzDaily/ ./ cd pitzDaily ./Allrun paraFoam
Cannot use ParaView reader module library (PVFoamReader) - not built?
Consider building the reader module
cd $WM_PROJECT_DIR/applications/utilities/postProcessing/graphics/PVReaders
./Allwclean
./Allwmake
Continuing with builtin reader: paraFoam -vtk
※Skip Zero Timeの表示を外すこと
エラーメッセージが表示されましたが、結果を表示することができました。
このエラーメッセージは「PVFoamReaderはインストールされませんでしたがvtkOpenFOAMReaderを使ってケースを表示しました」みたいな意味です。 PVFoamReaderはOpenFOAMで用意されたPluginで、vtkOpenFOAMReaderは元々ParaViewに備わっているVTK由来の機能です。 主な違いはSetを読み込めるかどうかかだったと思います。
なお今現在は表示されているエラーメッセージ通りにコマンドを入力してもPVFoamReaderのコンパイルは成功しません。（Qtの問題？）
SCOTCH(decomposePar)
次にThirdPartyとしてインストールされている並列計算のための分割手法であるSCOTCHの確認を行います。 以下のようにchtMultiRegionFoamのチュートリアルをコピーしてきて実行するとデフォルトでは並列数4で計算が行われます。 その際にdecomposeParではSCOTCHを用いて分割が行われます。
cd $WM_PROJECT_USER_DIR cp -r $FOAM_TUTORIALS/heatTransfer/chtMultiRegionFoam/multiRegionHeater/ ./ ./Allrun paraFoam -builtin
ちゃんと計算が行われていたらOKです。
Python3
最後にPython3の機能が搭載されているかをチェックします。
先程開いたParaViewのメニューバーでView→Python Shellをチェックすると下にコンソールが現れます。
ここをクリックすると以下のように環境が表示されますので、Python3環境であることが確認できます。
>>> Python 3.6.8 (default, Jan 14 2019, 11:02:34) [GCC 8.0.1 20180414 (experimental) [trunk revision 259383]] on linux >>> from paraview.simple import *
試しに以下のようなPython3のコマンドを実行してみても良いかと思います。
最後に
お疲れ様でした！
【OpenFOAM】ゴルフボールはどう打ったら最もよく飛ぶのか
はじめに
なかなかゴルフが上手くならないのでゴルフボールの周りの空気の流れを計算して最も飛距離を伸ばすための打ち出し角度と回転速度を求めてみました。
- バウンドの計算は着地点の角度や芝の状態などにより大きく変化すると考えられるため、着地するまでのキャリー距離を飛距離として計算していきます。
飛距離を追求してもゴルフの上達に繋がらない可能性があります。
2019.5.6 17:40修正 コメント頂いた回転行列の部分を修正しました。これにより最も飛距離が伸びるのは低角度でバックスピンという結果になりました。
方針
まず弾道を計算するにあたり以下の条件を固定します。
- ボールの初速は60m/s (ヘッドスピード45m/sで真芯に当たった場合)
- ボールの進行方向に対して左右方向への移動が無い
- ボールの回転方向もトップスピン/バックスピンの方向に固定
この状況で飛距離を計算するために以下のようにボールの弾道を計算する必要があります。
弾道計算を行う中では、空気抵抗とそれにより影響されるボールの速度や回転速度を同時に計算しつつ、その時のボールの位置を計算していくことになります。 この計算を行うには以下の方法が考えられます。
- 空気抵抗と弾道を同時に計算：
- 非定常ソルバーの中で、各時刻でボールにかかる力を求める↔境界条件を更新を繰り返す中で飛距離を計算するソルバーを作成する
- このソルバーによる飛距離計算を複数パターンで行い最も飛距離の長い条件を探索する。
- 空気抵抗と弾道を別々に計算：
- 定常ソルバーを用いて複数条件で定常状態のボールにかかる力を求める
- その結果から条件と力の関係を近似する
- その近似式による飛距離計算を複数パターンで行い最も飛距離の長い条件を探索する。
最終的に今回は2を採用しました。
今回のケースでは1の条件はあまりに計算に時間がかかり過ぎてしまい、PCのスペックの関係で断念しました。
今回のように速度に対して対象となるメッシュが小さい場合にはどうしても時間解像度が細かくなり計算負荷が大きくなってしまいます。
ただ本格的に計算する場合には非定常でのみ計算できる状況（カルマン渦など）の影響を無視することができないため、相応の環境を用意して1を計算する必要があるかと思います。
というわけで、具体的には
- OpenFOAMのsimpleFoamで回転するゴルフボール周りの定常流体解析
- そこで得られた抗力・揚力・モーメントを用いてPythonで弾道と飛距離を数値計算
という流れで飛距離を計算していきます。
計算環境
- OS : Ubuntu 18.04 LTS
- メッシャ：SALOME 8.4.0
- ソルバー：simpleFoam (OpenFOAM v1812)
- Python：3.6.6 (Anaconda 4.5.2)
OpenFOAMによる空気抵抗の計算
まずはOpenFOAMを用いてゴルフボール周りの流れを計算していきます。乱流にはk-ω SSTモデルを使用し、回転はMRFで計算しています。
メッシュ
ゴルフボールの規格を参考に以下のようなゴルフボールを想定しました。
- 直径：42.67 mm
- 重量：45.93 g
- ディンプル
- 数：362個
- 直径：1.8 mm
- 深さ：0.6 mm
SALOMEを用いて下図のようなメッシュを作成しました。
全てヘキサで作成した後にrefineMeshにて一部をポリヘドラに細分化しています。また水色と黄色のエリアは別々に作成しており、その境界はcyclicAMIにより接続されています。
$ checkMesh
Mesh stats
points: 2428023
faces: 6752074
cells: 2162148
Overall number of cells of each type:
hexahedra: 2108612
polyhedra: 53536
Checking geometry...
Overall domain bounding box (-0.2 -0.1 -0.1) (0.3 0.1 0.1)
Mesh has 3 geometric (non-empty/wedge) directions (1 1 1)
Mesh has 3 solution (non-empty) directions (1 1 1)
Boundary openness (-1.07072e-16 7.16432e-18 -1.11005e-15) OK.
Max cell openness = 3.42628e-16 OK.
Max aspect ratio = 4.95719 OK.
Minimum face area = 1.84412e-08. Maximum face area = 3.3335e-05. Face area magnitudes OK.
Min volume = 7.28758e-12. Max volume = 1.66675e-07. Total volume = 0.0199604. Cell volumes OK.
Mesh non-orthogonality Max: 55.035 average: 14.9621
Non-orthogonality check OK.
Face pyramids OK.
Max skewness = 0.987239 OK.
Coupled point location match (average 0) OK.
基本となるケースの作成
OpenFOAMのソルバーであるsimpleFoam用のケースディレクトリを作成します。
$FOAM_TUTORIALS/incompressible/simpleFoam/motorBikeとほぼ一緒の条件でメッシュとcontrolDictと緩和係数を変更しました。
ここでcontrolDict内ではfunctionObjectを使ってボールにかかる力を求めています。
/*--------------------------------*- C++ -*----------------------------------*\ | ========= | | | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox | | \\ / O peration | Version: v1812 | | \\ / A nd | Web: www.OpenFOAM.com | | \\/ M anipulation | | \*---------------------------------------------------------------------------*/ FoamFile { version 2.0; format ascii; class dictionary; object controlDict; } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // application simpleFoam; startFrom latestTime; startTime 0; stopAt endTime; endTime 20000; deltaT 1; writeControl timeStep; writeInterval 1000; purgeWrite 1; writeFormat binary; writePrecision 6; writeCompression off; timeFormat general; timePrecision 6; runTimeModifiable true; functions { force1 { type forces; libs ( "libforces.so" ); writeControl writeTime; log yes; patches ( wall_ball ); rho rhoInf; rhoInf 1; CofR (0 0 0); } } // ************************************************************************* //
また以下のように強めな緩和係数を設定しました。
SIMPLE
{
nNonOrthogonalCorrectors 1;
consistent yes;
}
potentialFlow
{
nNonOrthogonalCorrectors 10;
}
relaxationFactors
{
fields
{
p 0.4;
}
equations
{
U 0.4;
k 0.3;
omega 0.3;
}
}
条件毎にソルバー実行
今回は速度と回転軸を固定し、複数の回転速度で計算を行います。そのために以下のようにPythonにてケーススタディを行います。
まず以下のような関数を作成します。
import os, shutil from datetime import datetime import threading, time, glob, subprocess, requests import pandas as pd import paraview.simple as pv from PyFoam.RunDictionary.ParsedParameterFile import ParsedParameterFile from math import pi, sin, cos def makeGolfCase(resultDir, orgCase, **kwargs): if 'newName' in kwargs.keys(): baseName = kwargs['newName'] else: baseName = os.path.basename(orgCase) n = 0 newCase = os.path.join(resultDir, '{}_{}'.format(baseName, n)) while os.path.exists(newCase): newCase = os.path.join(resultDir, '{}_{}'.format(baseName, n)) n += 1 print('Created new case named', newCase) name = os.path.basename(newCase) shutil.copytree(orgCase,newCase) if not os.path.exists(newCase+'/0'): if os.path.exists(newCase+'/0.orig'): shutil.copytree(newCase+'/0.orig', newCase+'/0') if 'U' in kwargs: u = float(kwargs['U']) with ParsedParameterFile(newCase+'/0/U') as UFile: if kwargs['isInitialCase']: UFile.content['internalField'] = 'uniform ({:.2f} 0 0)'.format(u) UFile.content['boundaryField']['inlet']['value'] = 'uniform ({:.2f} 0 0)'.format(u) UFile.writeFile() with ParsedParameterFile(newCase+'/system/controlDict') as controlDict: controlDict.content['functions']['force1']['magUInf'] = '{:.2f}'.format(u) controlDict.writeFile() if 'rpm' in kwargs: if kwargs['rpm'] > 0: omega = kwargs['rpm'] * 2 * pi / 60.0 with ParsedParameterFile(newCase+'/constant/MRFProperties') as MRF: MRF.content['MRF1']['active'] = 'true' MRF.content['MRF1']['omega'] = 'constant {:.6f}'.format(omega) MRF.content['MRF1']['axis'] = '({:.6f} {:.6f} {:.6f})'.format(kwargs['axis'][0],kwargs['axis'][1],kwargs['axis'][2]) MRF.writeFile() return newCase
この中では
- ケースのコピー（名前が重複しないように）
- 速度Uを設定
- MRF回転条件を設定
を行っています。この関数を用いて以下のように実行します。
masterDir = '../study' if not os.path.exists(masterDir): os.mkdir(masterDir) # このkwargsをmakeGolfCase関数に渡す kwargs = { 'U' : 5, 'rpm' : 100, 'axis' : [0,0,1], 'isInitialCase' : True } for U in [40,50,60]: for rpm in [60, 1000, 3000, 5000 ]: for angle in [0]: kwargs['U'] = U kwargs['rpm'] = rpm kwargs['axis'] = [sin(angle*pi/180), 0, cos(angle*pi/180)] newName = makeGolfCase(masterDir, '../makeMesh/merge/case', **kwargs) cwdDir = os.getcwd() os.chdir(newName) subprocess.check_output(['./Allrun']) os.chdir(cwdDir)
これにより速度を[40,50,60]m/s × 回転速度[60, 1000, 3000, 5000]rpmの12パターンで計算を行います。
これで以下のようなボールにかかる力の計算が12パターン行われました。
Pythonによる弾道の計算
それでは得られた結果をもとに弾道と飛距離を計算していきます。使用するのはOpenFOAMのfunctionObjectsで出力されるファイルです。これをpandasで読み込んで加工しながら使っていきます。
結果の確認
計算が完了したらまずはちゃんと収束しているかどうかを確認してみましょう。 以下のようにfunctionObjectsにより生成されたforceとmomentをtimeStepを横軸にプロットしてみます。 今回は一度pandasに放り込むことで簡単に描画しています。
caseDir = '../study' cases = [] for d in glob.glob(caseDir+'/*'): if os.path.exists(d+'/system/controlDict'): cases.append(d) cases = sorted(cases) df = pd.DataFrame() fig = plt.figure(figsize=(20,10)) ind = 0 for case in cases: name = os.path.basename(case) forceFile = case+'/postProcessing/force1/0/force.dat' momdentFile = case+'/postProcessing/force1/0/moment.dat' forceDF = pd.read_csv(forceFile, sep='\t', header=3,index_col=0) force = np.array([ np.array(valStr[1:-1].split(' '), dtype=float) for key,valStr in forceDF[forceDF.columns[0]].items() ]) momentDF = pd.read_csv(momdentFile, sep='\t', header=3,index_col=0) moment = np.array([ np.array(valStr[1:-1].split(' '), dtype=float) for key,valStr in momentDF[momentDF.columns[0]].items() ]) for i in range(3): plt.subplot(2,3,i+1) plt.plot(range(1000),force[-1000:,i], label=name) plt.subplot(2,3,i+4) plt.plot(range(1000),moment[-1000:,i], label=name) i = 0 for axis in ['_X','_Y','_Z']: plt.subplot(2,3,i+1) plt.title('force'+axis) plt.subplot(2,3,i+4) plt.title('moment'+axis) if i == 2: plt.legend() i += 1 fig.savefig('fig/1.results.png', bbox_inches="tight")
force_Xとforce_Yとmoment_Zはいい感じに収束しています。 なお今回の計算条件ではZ軸が回転軸となっているので
- X方向：抗力
- Y方向：揚力
- Z軸モーメント：回転モーメント
となります。
そのため今回の弾道計算の中では進行方向に対して左右にかかる力(force_Z)は無視します。
またモーメントについても回転を弱める方向のモーメント(moment_Z)以外については無視することとします。
というわけでよく収束されていたforce_Xとforce_Yとmoment_Zのみを用いて弾道を計算していきます。
各条件における計算結果を抽出
各条件によりかかる力/モーメントを同じくpandasのDataFrame形式にまとめていきます。 今回は定常計算の各timeStepで得られる数値の後半100step分の結果を平均した値を結果として使用します。
df = pd.DataFrame() ind = 0 for case in cases: name = os.path.basename(case) forceFile = case+'/postProcessing/forceCoeffs1/0/force.dat' momdentFile = case+'/postProcessing/forceCoeffs1/0/moment.dat' forceDF = pd.read_csv(forceFile, sep='\t', header=3,index_col=0) force = np.array([ np.array(valStr[1:-1].split(' '), dtype=float) for key,valStr in forceDF[forceDF.columns[0]].items() ]) momentDF = pd.read_csv(momdentFile, sep='\t', header=3,index_col=0) moment = np.array([ np.array(valStr[1:-1].split(' '), dtype=float) for key,valStr in momentDF[momentDF.columns[0]].items() ]) rpm = int(round(ParsedParameterFile(case+'/constant/MRFProperties').content['MRF1']['omega'][1] * 60 / 2 / pi, 0)) U = int(round(ParsedParameterFile(case+'/system/controlDict').content['functions']['forceCoeffs1']['magUInf'])) ser = {'rpm':rpm,'U':U} n = 0 for axis in ['_X','_Y','_Z']: ser['force'+axis] = force[-100:,n].mean() n += 1 n = 0 for axis in ['_X','_Y','_Z']: ser['moment'+axis] = moment[-100:,n].mean() n += 1 df[ind] = pd.Series(ser) ind += 1 # 逆回転 ser = {'rpm':-rpm,'U':U} n = 0 for axis in ['_X','_Y','_Z']: if n != 1: ser['force'+axis] = force[-int(len(force)/10):,n].mean() else: ser['force'+axis] = -force[-int(len(force)/10):,n].mean() n += 1 n = 0 for axis in ['_X','_Y','_Z']: if n != 2: ser['moment'+axis] = moment[-int(len(moment)/10):,n].mean() else: ser['moment'+axis] = -moment[-int(len(moment)/10):,n].mean() n += 1 df[ind] = pd.Series(ser) ind += 1 df = df.T df = df.set_index(['rpm','U'], drop=True) df
定常計算の中では重力加速度については計算されていないため、Y方向を反転させた結果をそのまま逆回転の結果として使用しています。 あえてこの結果をDataFrameに入れた理由は、後に出てくる近似式の無回転状態の精度を高めるためです。
データの傾向を確認
以上のようにまとめた結果を以下のように図示してみます。
fig = plt.figure(figsize=(20,10)) ind = 0 for name in ['force_X', 'force_Y', 'moment_Z']: ans = {'rpm':[],'U':[]} for key,val in df[name].items(): rpm = key[0] U = key[1] rpmKey = '{}rpm'.format(int(rpm)) UKey = '{}m/s'.format(int(U)) if not rpmKey in ans: ans[rpmKey] = [] ans['rpm'].append(rpm) if not UKey in ans: ans[UKey] = [] ans['U'].append(U) ans[rpmKey].append(val) ans[UKey].append(val) plt.subplot(2,3,ind+1) for rpm in ans['rpm']: rpmKey = '{}rpm'.format(int(rpm)) plt.scatter(ans['U'],ans[rpmKey],label=rpmKey) plt.title(name) plt.legend() plt.subplot(2,3,ind+4) for U in ans['U']: UKey = '{}m/s'.format(int(U)) plt.scatter(ans['rpm'],ans[UKey],label=UKey) plt.title(name) plt.legend() ind += 1 fig.savefig('fig/2.plot.png', bbox_inches="tight")
なんとなく綺麗な規則性がありそうですね。 少々雑ですが以下のような式で近似してみましょう。
- U : 速度 [m/s]
- rot : 回転速度 [rpm]
- a,b,c,d,e,f : 定数
近似にはscipyのcurve_fitを用います。ます以下のような関数を用意します。
from scipy.optimize import curve_fit def func(xy, a,b,c,d,e,f): x = xy[0] y = xy[1] return (a*x**2 + b*x + c)*(d*y**2+e*y+f)
この上で以下のように実行することでforce_X、force_Y、moment_Zそれぞれについてフィッティングが行われます。
fig = plt.figure(figsize=(20,20)) rpmX = list(range(-5000,5001,100)) UX = list(range(0,61,5)) betas = {} ind = 0 for name in ['force_X', 'force_Y', 'moment_Z']: ans = {'rpm':[0.0],'U':[0.0]} for key,val in df[name].items(): rpm = key[0] U = key[1] rpmKey = '{}rpm'.format(int(rpm)) UKey = '{}m/s'.format(int(U)) if not rpmKey in ans: ans[rpmKey] = [0.0] ans['rpm'].append(rpm) if not UKey in ans: ans[UKey] = [0.0] ans['U'].append(U) ans[rpmKey].append(val) ans[UKey].append(val) ans['0rpm'] = [v for v in ans[rpmKey]] ans['0m/s'] = [v for v in ans[UKey]] beta, pcov = curve_fit(func, np.array(list(df[name].index)).T, df[name].values) betas[name] = beta plt.subplot(3,3,ind+1) for rpm in ans['rpm']: rpmKey = '{}rpm'.format(int(rpm)) plt.plot(UX, [func((rpm,u),beta[0],beta[1],beta[2],beta[3],beta[4],beta[5]) for u in UX], linestyle='dashed',label='fitting curve') plt.scatter(ans['U'],ans[rpmKey],label=rpmKey) plt.title(name+' vs velocity') plt.xlabel('velocity [m/s]') plt.legend() plt.subplot(3,3,ind+4) for U in ans['U']: UKey = '{}m/s'.format(int(U)) plt.plot(rpmX, [func((r,U),beta[0],beta[1],beta[2],beta[3],beta[4],beta[5]) for r in rpmX], linestyle='dashed',label='fitting curve') plt.scatter(ans['rpm'],ans[UKey],label=UKey) plt.title(name+' vs rotation') plt.xlabel('rotation [rpm]') plt.legend() rpms, Us = np.mgrid[-5000:5000:100, 0:61:5] Z = func((rpms, Us),beta[0],beta[1],beta[2],beta[3],beta[4],beta[5]) plt.subplot(3,3,ind+7) plt.contour(rpms, Us, Z, linewidths=0.5, colors='k') plt.contourf(rpms, Us, Z, cmap=plt.cm.Spectral) plt.title(name+' contour') plt.xlabel('rotation [rpm]') plt.ylabel('velocity [m/s]') ind += 1 fig.savefig('fig/3.curvefit.png', bbox_inches="tight")
この中の下のコンター図は近似式を用いて計算した値をプロットしたものです。 OpenFOAMで計算していない条件についてもいい感じに補完されているように見えます。
なお得られた近似式の各定数は以下のようになりました。
tmpdf = pd.DataFrame(betas).T tmpdf.columns = ['a', 'b', 'c', 'd', 'e', 'f'] tmpdf
以上のようにOpenFOAMの結果に則った近似式がそれぞれの力についてできました。
近似式を用いて弾道を計算する
速度と回転速度を与えたらボールに加わる力が得られるようになりましたのでこれを使って弾道を計算していきます。
弾道計算にはオイラー法っぽい感じで時刻を細かく刻んで進行していき
前の時刻の結果から次の時刻の条件を計算するといったことを繰り返していきます。
2019.5.6 17:40修正 コメントをいただき回転行列の箇所を修正しました。
r = 42.67e-3/2 # m m = 45.93e-3 # kg def newStatus(rpms, vecU, dt, pos): pos.append([pos[-1][0]+vecU[0]*dt, pos[-1][1]+vecU[1]*dt]) rpm = rpms[-1] magU = (vecU[0]**2 + vecU[1]**2)**0.5 cosT = vecU[0] / magU sinT = vecU[1] / magU f = {} for key,beta in betas.items(): f[key] = func((rpm, magU),beta[0],beta[1],beta[2],beta[3],beta[4],beta[5]) force = [-f['force_X']*cosT - f['force_Y']*sinT, -f['force_X']*sinT + f['force_Y']*cosT] # 修正しました。 #force = [-f['force_X']*sinT - f['force_Y']*cosT, -f['force_X']*cosT + f['force_Y']*sinT] vecU = [vecU[0] + force[0]/m*dt, vecU[1] + (force[1]/m - 9.81)*dt] rpms.append(rpm + f['moment_Z']/((2/5)*m*r**2)*dt*60/2/pi) #print(vecU, pos[-1]) if pos[-1][1] > 0 and len(pos) < 1000: newStatus(rpms, vecU, dt, pos) def distance(magU, rpm, angle): pos = [[0.0, 0.0]] vecU = [magU*cos(pi*angle/180),magU*sin(pi*angle/180)] rpms =[rpm] dt = 1e-3 newStatus(rpms, vecU, dt, pos) return np.array(pos)[:,0].max() def distancePlot(magU, rpm, angle): pos = [[0.0, 0.0]] vecU = [magU*cos(pi*angle/180),magU*sin(pi*angle/180)] rpms =[rpm] dt = 1e-3 topback = 'top' if rpm > 0: topback = 'back' label = '{}m/s, {}deg, {}rpm {} spin'.format(int(magU), int(angle), int(abs(rpm)), topback) newStatus(rpms, vecU, dt, pos) p = np.array(pos) plt.subplot(2,1,1) plt.plot(p[:,0],p[:,1],label=label) plt.ylabel('height [m]');plt.xlim(0,200);plt.ylim(0,100) plt.subplot(2,1,2) plt.plot(p[:,0],rpms,label=label) plt.xlabel('length [m]');plt.ylabel('rotating speed [rpm]'); plt.xlim(0,200);plt.ylim(-5000,5000)
説明の前に実行してみましょう。以下のように速度一定=60m/s、仰角一定=45°、回転速度を[-5000,-2500,0,2500,5000]の5パターンで計算を行うと各弾道を得ることができます。
fig = plt.figure(figsize=(15,7)) for rpm in [-5000,-2500,0,2500,5000]: distancePlot(60, rpm, 45) plt.legend() fig.savefig('fig/4.lines.png', bbox_inches="tight")
2019.5.6 17:40修正 図を差し替えました。
この上の図が弾道です。横軸に距離、縦軸に高さを取っています。また下の図では回転速度が減衰していっている様子が図示されています。
newStatus関数は再帰関数になっており、posというリストに各時刻の位置がappendされていきます。 この各時刻において速度と回転速度から力が計算され、そこからさっきの近似式を用いて次の時刻の回転速度と速度ベクトルが求められます。
同じような内容のdistance関数では飛距離を返すようになっています。
複数条件で計算
以上のようにdistance関数で各条件で飛距離を求めることができるようになりましたので、多くのパターンで計算を行ってみましょう。 なお速度は60m/sに固定した条件で計算を行います。 また最終的に2変数のコンター図にしたいのでX(回転速度)、Y(仰角)、Z(飛距離)で結果をまとめられるように計算を行います。
X, Y = np.mgrid[0:8001:200, 6:61:2] Z = np.zeros_like(X) maxDist = 0.0 maxCond = [0,0] for i in range(len(X)): for j in range(len(X[i])): val = distance(60, X[i][j], Y[i][j]) Z[i][j] = val if val > maxDist: maxDist = val maxCond = [X[i][j], Y[i][j]] print('Maximum distance is {:.1f} m when {:.1f} rpm and {:.1f} degree'.format(maxDist, maxCond[0], maxCond[1])) fig = plt.figure(figsize=(15,10)) plt.contour(X, Y, Z, linewidths=0.5, colors='k') plt.pcolor(X,Y, Z, cmap=plt.cm.Spectral) plt.xlabel('rotation [rpm] (+ : back spin, - : top spin)') plt.ylabel('shot angle [degree]') plt.colorbar() fig.savefig('fig/5.pcolor.png', bbox_inches="tight")
2019.5.6 17:40修正 図を差し替えました。
Maximum distance is 180.2 m when 6000.0 rpm and 18.0 degree
という訳で60m/sで打った時に最も飛距離がでるのは仰角18°の方向に6000rpmのバックスピンであることがわかりました。
最も飛ぶ弾道を図示
最後にこの条件を図示してみましょう。
fig = plt.figure(figsize=(15,7)) distancePlot(60, maxCond[0], maxCond[1]) plt.legend() fig.savefig('fig/6.mostFar.png', bbox_inches="tight")
2019.5.6 17:40修正 図を差し替えました。
低めに高回転でホップさせる感じがいいのでしょうか。
最後に
いかがでしたか？
みなさんもドラコンなどの際には仰角18°の方向に6000rpmのバックスピンをかけて打ってみてください。
【SALOME】スクリプトでメッシュのケーススタディ
はじめに
以前にSALOMEを使ったヘキサメッシュの作成方法を紹介しました。
これはSALOMEのGUIを用いた方法ですが、SALOMEはTUIの中で使用すると更に便利だったりします。 今回はTUIでちょっと寸法を変えたメッシュを量産していく方法を紹介していきます。
お試しスクリプト
例として以下のようなスクリプトを用意します。
sampleBox.py
# -*- coding: utf-8 -*- import os, sys, salome, GEOM, SMESH from salome.geom import geomBuilder from salome.smesh import smeshBuilder # スタディとGEOM(形状)およびSMESH(メッシュ)のオブジェクトを作成 salome.salome_init() theStudy = salome.myStudy geompy = geomBuilder.New(theStudy) smesh = smeshBuilder.New(theStudy) # 1m×1m×1mの箱を作成 Box_1 = geompy.MakeBoxDXDYDZ(1, 1, 1) geompy.addToStudy( Box_1, 'Box_1' ) # 各面にグループ名をつける groupNo = {'minZ':31, 'maxZ':33, 'minY':23, 'maxY':27, 'minX':3, 'maxX':13} geomGroups = {} for key, i in groupNo.items(): geomGroups[key] = geompy.CreateGroup(Box_1, geompy.ShapeType["FACE"]) geompy.UnionIDs(geomGroups[key], [i]) geompy.addToStudyInFather( Box_1, geomGroups[key], key ) # メッシュを作成 Mesh_1 = smesh.Mesh(Box_1) Regular_1D = Mesh_1.Segment() Number_of_Segments_1 = Regular_1D.NumberOfSegments(15) Quadrangle_2D = Mesh_1.Quadrangle(algo=smeshBuilder.QUADRANGLE) Hexa_3D = Mesh_1.Hexahedron(algo=smeshBuilder.Hexa) isDone = Mesh_1.Compute() # Create groups from geometry for key, group in geomGroups.items(): tmp = Mesh_1.GroupOnGeom(group,key,SMESH.FACE) # UNVファイルとしてエクスポート Mesh_1.ExportUNV( 'Mesh.unv' )
このスクリプトは1m×1m×1mの箱を作成し各辺を15分割した直方体メッシュを生成するものです。 SALOMEのGUIでFile→Load Scriptでこのファイルを読み込むと以下のようなメッシュがGUIでできていることを確認できるかと思います。 また$HOMEに"Mesh.unv"が生成されているかと思います。
SALOMEのTUIモード
このスクリプトをSALOMEを起動せずに実行するにはSALOMEをインストールしたディレクトリにある"salome"ファイルに引数" -w 1 -t "を付けて以下のようにターミナルに入力します。
$ path/to/salome -w 1 -t sampleBox.py
runSalome running on inabower-LAPTOP
Searching for a free port for naming service: 2812 - OK
Searching Naming Service + found in 0.1 seconds
Searching /Registry in Naming Service + found in 0.5 seconds
Searching /Kernel/ModulCatalog in Naming Service +th. 140452740915712 - Trace /home/rd-ap-palmco/salome/edf/V8_4_BR/V8_4_0/modules/src/KERNEL/src/ModuleCatalog/SALOME_ModuleCatalog_Server.cxx [101] : Module Catalog Server: Naming Service was found
Warning: this type (Study,objref) already exists, it will be ignored.
Warning: this type (pyobj,objref) already exists, it will be ignored.
Warning: this type (SALOME_MED/MEDCouplingFieldDoubleCorbaInterface,objref) already exists, it will be ignored.
found in 0.5 seconds
RunStudy
Searching /myStudyManager in Naming Service + found in 0.5 seconds
Searching /Containers/inabower-LAPTOP/FactoryServer in Naming Service + found in 0.5 seconds
Start SALOME, elapsed time : 2.1 seconds
createNewStudy
extStudy 1
※なおAnacondaなどでデフォルトをPython3にしている場合には、Python2で書かれたSALOMEを実行することができないためPython2の環境をactivateしておく必要があります。
これで今度はsampleScript.pyと同じ場所に"Mesh.unv"が生成されていることが確認できるかと思います。
関数によるコマンドの省略
ただ毎回SALOMEの場所を入力するのは手間ですので~/.bashrcに関数やaliasを作成しておくと良いかと思います。
~/.bashrc
function runSalome (){ path/to/salome -w 1 -t $1 $2 } export -f runSalome
Anacondaで毎回Python2にするのがしんどいという方は以下のようにすると良いかと思います。
~/.bashrc
function runSalome (){ source activate py2 # 各自のpython2環境 path/to/salome -w 1 -t $1 $2 source activate base # 元の環境に戻す } export -f runSalome
これで以下のコマンドで同様の実行を行うことができるようになりました。
$ runSalome sampleBox.py
runSalome running on inabower-LAPTOP
Searching for a free port for naming service: 2812 - OK
Searching Naming Service + found in 0.1 seconds
Searching /Registry in Naming Service + found in 0.5 seconds
Searching /Kernel/ModulCatalog in Naming Service +th. 140452740915712 - Trace /home/rd-ap-palmco/salome/edf/V8_4_BR/V8_4_0/modules/src/KERNEL/src/ModuleCatalog/SALOME_ModuleCatalog_Server.cxx [101] : Module Catalog Server: Naming Service was found
Warning: this type (Study,objref) already exists, it will be ignored.
Warning: this type (pyobj,objref) already exists, it will be ignored.
Warning: this type (SALOME_MED/MEDCouplingFieldDoubleCorbaInterface,objref) already exists, it will be ignored.
found in 0.5 seconds
RunStudy
Searching /myStudyManager in Naming Service + found in 0.5 seconds
Searching /Containers/inabower-LAPTOP/FactoryServer in Naming Service + found in 0.5 seconds
Start SALOME, elapsed time : 2.1 seconds
createNewStudy
extStudy 1
引数の追加
このスクリプトには引数をつけることもできます。 ヘルプを見てみましょう。
$ runSalome -help
Usage: salome start [options] [STUDY_FILE] [PYTHON_FILE [args] [PYTHON_FILE [args]...]]
Python file arguments, if any, must be comma-separated (without blank characters) and prefixed by "args:" (without quotes), e.g. myscript.py args:arg1,arg2=val,...
Options:
-t, --terminal Launch without GUI (in the terminal mode).
-g, --gui Launch in GUI mode [default].
(以下略)
というわけで"args:aaa,bbb"のように引数を入力すれば良いことが分かります。 引数で箱の大きさと出力するファイル名を入力できるようにしてみます。
sampleBoxArgs.py
# -*- coding: utf-8 -*- import os, sys, salome, GEOM, SMESH from salome.geom import geomBuilder from salome.smesh import smeshBuilder # スタディとGEOM(形状)およびSMESH(メッシュ)のオブジェクトを作成 salome.salome_init() theStudy = salome.myStudy geompy = geomBuilder.New(theStudy) smesh = smeshBuilder.New(theStudy) # 引数を読み込む print sys.argv try: length = float(sys.argv[1]) fileName = sys.argv[2] if not fileName.endswith('.unv'): raise Exception('File name must be *.unv file') except: raise Exception('Wrong args') # 1m×1m×1mの箱を作成 Box_1 = geompy.MakeBoxDXDYDZ(length,length,length) geompy.addToStudy( Box_1, 'Box_1' ) # 各面にグループ名をつける groupNo = {'minZ':[31], 'maxZ':[33], 'minY':[23], 'maxY':[27], 'minX':[3], 'maxX':[13]} geomGroups = {} for key, faces in groupNo.items(): geomGroups[key] = geompy.CreateGroup(Box_1, geompy.ShapeType["FACE"]) geompy.UnionIDs(geomGroups[key], faces) geompy.addToStudyInFather( Box_1, geomGroups[key], key ) # メッシュを作成 Mesh_1 = smesh.Mesh(Box_1) Regular_1D = Mesh_1.Segment() Number_of_Segments_1 = Regular_1D.NumberOfSegments(15) Quadrangle_2D = Mesh_1.Quadrangle(algo=smeshBuilder.QUADRANGLE) Hexa_3D = Mesh_1.Hexahedron(algo=smeshBuilder.Hexa) isDone = Mesh_1.Compute() # Create groups from geometry for key, group in geomGroups.items(): tmp = Mesh_1.GroupOnGeom(group,key,SMESH.FACE) # UNVファイルとしてエクスポート Mesh_1.ExportUNV( fileName )
これで以下のように入力すると各辺2mの箱のメッシュが"test.unv"という名前で生成されます。
$ runSalome sampleBoxArgs.py args:2,test.unv
runSalome running on inabower-LAPTOP
Searching for a free port for naming service: 2818 - OK
Searching Naming Service + found in 0.1 seconds
Searching /Registry in Naming Service + found in 0.5 seconds
Searching /Kernel/ModulCatalog in Naming Service +th. 140346892026368 - Trace /home/rd-ap-palmco/salome/edf/V8_4_BR/V8_4_0/modules/src/KERNEL/src/ModuleCatalog/SALOME_ModuleCatalog_Server.cxx [101] : Module Catalog Server: Naming Service was found
Warning: this type (Study,objref) already exists, it will be ignored.
Warning: this type (pyobj,objref) already exists, it will be ignored.
Warning: this type (SALOME_MED/MEDCouplingFieldDoubleCorbaInterface,objref) already exists, it will be ignored.
found in 0.5 seconds
RunStudy
Searching /myStudyManager in Naming Service + found in 0.5 seconds
Searching /Containers/inabower-LAPTOP/FactoryServer in Naming Service + found in 0.5 seconds
Start SALOME, elapsed time : 2.2 seconds
createNewStudy
extStudy 1
['/home/inabower/Documents/Models/sampleBox/sampleBoxArgs.py', '2', 'test.unv']
というわけでTUIでスクリプトを実行できるようになりました。
ケーススタディ
それでは辺の長さを1m,2m,3m,4m,5mで一括メッシュ生成をしてみます。 以下のようなシェルスクリプトを作ります。
#!/bin/bash mkdir mesh for n in 1 2 3 4 5; do echo "Making mesh. Length = ${n} m" runSalome sampleBoxArgs.py "args:${n},mesh/test${n}.unv" done
これで以下のようにメッシュが一気に生成されました。
$ tree mesh
mesh
├── test1.unv
├── test2.unv
├── test3.unv
├── test4.unv
└── test5.unv
0 directories, 5 files
なおOpenFOAMのケースを一気に作る場合は、まずorgという元となるケースを準備します。
tree org
org
├── makeMesh.sh
├── constant
└── system
├── controlDict
├── fvSchemes
└── fvSolution
2 directories, 4 files
makeMesh.shの中身はこんな感じ。
#!/bin/bash
ideasUnvToFoam box.unv
その上で以下のようなスクリプトでunvの変換までを一気にやります。
#!/bin/bash mkdir mesh for n in 1 2 3 4 5; do echo "Making mesh. Length = ${n} m" cp -r org mesh/test${n} runSalome sampleBoxArgs.py "args:${n},mesh/test${n}/box.unv" pushd mesh/test${n} ./makeMesh.sh touch test${n}.foam popd done
さらにソルバーやpvbatchによる図の出力などを一括でできるようになるととても便利です。
スクリプトの作り方
SALOMEでメッシュを作成した後にFile→Dump Studyでそれまで行った作業をスクリプトにしてくれる機能があります。 最初は「GUIで少し作業→スクリプトに吐き出す→中身を少し改造する→GUIで読み込む」を繰り返すことでスクリプトの作り方がわかってくるかと思います。 この辺りもそのうち掘り下げて説明できたらと思います。
GUIでの作業方法は以下の記事で説明しています。
またこんな機能無いかなと思ったら以下の公式APIを見てみましょう。
最後に
SALOMEはとても自由度の高いメッシャですので是非試してみてください。
