こんにちは（@t_kun_kamakiri）
前回の記事では粒子追跡の設定について解説しました。
今回は流速と粒子直径をパラメータとしたパラメータスタディをPythonを用いて実行する方法を解説します。
Pythonを用いて粒子追跡解析（OpenFOAM）のパラメータスタディを行う
この記事を読むことでPythonを使ったOpenFOAMのパラメータスタディの基礎がわかるでしょう。
本記事のケースフォルダはgithubに公開しています。
OpenFOAM v2412(WSL2)
解析の対象
解析の対象はキャビティ流れで作った流れに乗って粒子が運動する粒子追跡解析モデルです。
キャビティ流れは流速を求めるための非定常解析を行い、得られた最終時刻の流速に乗って粒子が運動するというone-wayカップリングです。
パラスタのフォルダ構成
最終的なフォルダ構成は以下となっています。
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 |
tree -L 3 . ├── orgCase │ ├── 001_cavity │ │ ├── 0 │ │ ├── Allrun.icoFoam │ │ ├── constant │ │ └── system │ ├── 002_particle │ │ ├── 0 │ │ ├── Allrun.icoUncoupledKinematicParcelFoam │ │ ├── constant │ │ ├── post.foam │ │ └── system │ ├── Allrun │ └── orgCase.foam ├── parastduty.ipynb |
- orgCase：オリジナルのフォルダ
- 001_cavity：キャビティ流れ
- 002_particle：粒子追跡
このようにオリジナルのフォルダにはキャビティ流れ解析用と粒子追跡解析用の2つのフォルダを用意して別々で解析する構成にしています。
オリジナルモデルの作成
パラメータスタディするにあたっては、オリジナルのフォルダを計算できる状態に整えておく必要があります。
詳細は以下の記事に書いておりますが、本記事でも復習として簡単に流れを書いておきます。
パラメータスタディ用のフォルダを作成し、フォルダ移動を行います。
|
1 2 |
mkdir 20250209_paraStudy_particle cd 20250209_paraStudy_particle |
オリジナルフォルダを作成し、フォルダ移動を行います。
|
1 2 |
mkdir orgCase cd orgCase/ |
次にキャビティ流れと粒子追跡用のチュートリアルをコピーします。
|
1 2 |
cp -r $FOAM_TUTORIALS/incompressible/icoFoam/cavity/cavity 001_cavity cp -r /$FOAM_TUTORIALS/lagrangian/icoUncoupledKinematicParcelFoam/hopper/hopperInitialState 002_particle |
キャビティ流れの計算
まずはキャビティ流れの計算実行ができる状態にします。
フォルダ移動を行いblockMeshでメッシュ作成だけを行います。
|
1 2 |
cd 001_cavity blockMesh |
キャビティ流れのAllrun.icoFoamスクリプトを作成します。
Allrun.icoFoam
|
1 2 3 |
#!/bin/bash icoFoam |
これでキャビティ流れができる状態になります。
粒子追跡の計算
粒子追跡の方にフォルダ移動します。
|
1 |
cd ../002_particle |
メッシュ情報が異なるのでキャビティ流れのメッシュ情報をコピーします。
|
1 |
cp -r ../001_cavity/constant/polyMesh constant/ |
またキャビティ流れの境界条件もコピーします。
|
1 |
cp -r ../001_cavity/0 . |
キャビティ流れの0.5秒での流速分布と圧力分布をマッピングします。
|
1 |
mapFields ../001_cavity -sourceTime 0.5 -consistent |
mapFieldsによってキャビティ流れの最終状態を今から解析する粒子追跡の流速分布の初期状態として使うことができます。
ParaViewで初期状態の流速分布を見てみましょう。
キャビティ流れの流速分布が初期状態としてマッピングされています。
次に粒子の配置について設定を行います。
constant/kinematicCloudPositions
|
1 2 3 4 5 |
( (0.025 0.08 0.005) (0.050 0.08 0.005) (0.075 0.08 0.005) ) |
続いて粒子の設定についてです。
ここでは、粒子の密度や力の種類、壁との接触の定義などを行います。細かい設定内容についてはこだわらずほとんどデフォルトのままいきます。
constant/kinematicCloudProperties
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 |
solution { active true; coupled true; transient yes; cellValueSourceCorrection off; interpolationSchemes { rho cell; U cellPoint; mu cell; } integrationSchemes { U Euler; } sourceTerms { /* //U のスキームには "explicit" と "semiImplicit" が選べる。 スキームの後の数値は緩和係数らしい。 */ schemes { U explicit 1; } } } constantProperties { parcelTypeId 1; rhoMin 1e-15; minParticleMass 1e-15; rho0 964; youngsModulus 6e8; poissonsRatio 0.35; constantVolume false; } subModels { particleForces { sphereDrag; gravity; } injectionModels { model1 { type manualInjection; massTotal 0; parcelBasisType fixed; nParticle 1; SOI 0; positionsFile "kinematicCloudPositions"; U0 ( 0 0 0 ); sizeDistribution { type fixedValue; fixedValueDistribution { value 0.02; } } } } dispersionModel none; patchInteractionModel standardWallInteraction; heatTransferModel none; surfaceFilmModel none; stochasticCollisionModel none; collisionModel pairCollision; radiation off; pairCollisionCoeffs { // Maximum possible particle diameter expected at any time maxInteractionDistance 0.02; writeReferredParticleCloud no; pairModel pairSpringSliderDashpot; pairSpringSliderDashpotCoeffs { useEquivalentSize no; alpha 0.12; b 1.5; mu 0.52; cohesionEnergyDensity 0; collisionResolutionSteps 12; }; wallModel wallLocalSpringSliderDashpot; wallLocalSpringSliderDashpotCoeffs { useEquivalentSize no; collisionResolutionSteps 12; movingWall { youngsModulus 1e10; poissonsRatio 0.23; alpha 0.12; b 1.5; mu 0.43; cohesionEnergyDensity 0; } fixedWalls { youngsModulus 1e10; poissonsRatio 0.23; alpha 0.12; b 1.5; mu 0.43; cohesionEnergyDensity 0; } frontAndBack { youngsModulus 1e10; poissonsRatio 0.23; alpha 0.12; b 1.5; mu 0.1; cohesionEnergyDensity 0; } }; } standardWallInteractionCoeffs { type rebound; } } cloudFunctions {} |
設定が終わりましたので、計算実行します。
※必要に応じて「system/controlDict」で計算時間の設定を変えれば良いです。
Allrun.icoUncoupledKinematicParcelFoamのスクリプトを以下とします。
Allrun.icoUncoupledKinematicParcelFoam
|
1 2 3 4 5 |
#!/bin/bash mapFields ../001_cavity -sourceTime 0.5 -consistent icoUncoupledKinematicParcelFoam foamToVTK |
以上でキャビティ流れの流速に乗って運動する粒子追跡の計算が実行できるようになりました。
Allrunスクリプト
では、001_cavityと002_particleを連続で実行できるようにします。
orgCaseフォルダに移動します。
|
1 |
cd ../ |
そして、Allrunスクリプトを以下とします。
|
1 2 3 4 5 6 7 8 9 10 11 12 |
``` #!/bin/bash # cavityの計算 cd 001_cavity ./Allrun.icoFoam # 粒子追跡の計算 cd .. cd 002_particle ./Allrun.icoUncoupledKinematicParcelFoam ``` |
試しにAllrunを実行して正常に動作するかを確かめてください。
ParaViewは確認すると以下のようになっているでしょう。
PyFoamによるパラメータスタディ
ここからPythonを使ったOpenFOAMのパラメータスタディを行います。
全体のプログラムは以下としています。
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 |
import os import shutil import threading from pathlib import Path import subprocess from PyFoam.RunDictionary.SolutionDirectory import SolutionDirectory from PyFoam.RunDictionary.ParsedParameterFile import ParsedParameterFile from PyFoam.Execution.BasicRunner import BasicRunner def clone_file(orgCase, resultDir): "orgCaseをresultDirに重複しないようにコピー" baseName = Path(orgCase).name addName = None if addName!=None: baseName += f'_{addName}' Path(resultDir).mkdir(exist_ok=True, parents=True) n=0 newCase = Path(resultDir) / f'{baseName}_{n}' while newCase.is_dir(): n += 1 newCase = Path(resultDir) / f'{baseName}_{n}' shutil.copytree(orgCase, newCase) return newCase orgCase = "orgCase" resultDir = "results" U_list = [0.5, 1.0, 2.0, 3.0, 4.0, 5.0] d_mean_list = [0.005, 0.01, 0.02, 0.03, 0.04, 0.05] for U in U_list: for d_mean in d_mean_list: # 新しいファイルをコピー newCase = clone_file(orgCase, resultDir) # 流速の変更(001_cavity) cavity_dir = f"{newCase}/001_cavity" UFile = ParsedParameterFile(os.path.join(cavity_dir, "0", "U")) UFile["boundaryField"]["movingWall"]["value"] = f'uniform ({U} 0 0)' # 粒子径の変更 (002_particle) particle_dir = f"{newCase}/002_particle" kinematicCloudPropertiesFile = ParsedParameterFile(os.path.join(particle_dir, "constant", "kinematicCloudProperties")) kinematicCloudPropertiesFile["subModels"]["injectionModels"]["model1"]["sizeDistribution"]["fixedValueDistribution"]["value"] = d_mean kinematicCloudPropertiesFile.writeFile() print(f"{newCase}：全てのクローン操作とファイル編集が完了しました。") with open(os.devnull, 'w') as devnull: subprocess.run([f'./Allrun'], cwd=newCase, stdout=devnull, stderr=devnull) |
ここでは、キャビティ流れの流速U_list = [0.5, 1.0, 2.0, 3.0, 4.0, 5.0]と粒子直径 d_mean_list = [0.005, 0.01, 0.02, 0.03, 0.04, 0.05]をパラメータにしています。
本記事ではjupyter labを起動して実行しています。
技術書典でパラメータスタディ用の解説あり
技術書典17で出典でもOpenFOAMのパラスタの解説を行っていますので、是非手に取ってみてください。
内容はブログ書いているOpenFOAMの内容を元に、対象を「2次元円柱内流れ」として、
- カルマン渦の生成
- レイノルズ数と抗力係数の関係
- レイノルズ数とストローハル数
- ポテンシャル流れと粘性流体の圧力係数の違い
- 力学的相似性
- Pythonを使ったパラメータスタディ
などなど、モデル作成からグラフ作成までの一連の流れをまとめたものになります。
ページ数は140ページです。
↓内容のチラ見せ
まとめ
今回の記事では、OpenFOAM を用いた粒子追跡解析において、流速と粒子直径という2つのパラメータを変化させたパラメータスタディの実施方法を、PythonとPyFoamを使って自動化する手法を解説しました。
本記事のケースフォルダはgithubに公開しています。
次回の記事では、PyVistaを用いて、生成された VTK ファイルから解析結果の画像を自動生成するプログラムについて解説する予定です。
おすすめ参考図書
☟こちらは、OpenFOAMの日本語書籍が無い中唯一わかりやすい参考書だと思います。
☟こちらもOpenFOMの古いバージョンでの和訳になります。さすがにこちらはバージョンに対する日本語でのケアはしていないので、OpenFOAMに慣れている方は購入しても良いかと思います。僕は「日本語でまとまっている内容」なので少し重宝しています。
☟以下に、もっと初心者向けの同人誌を紹介しておきます。
初心者は「ってか、まずどうやってOpenFOAMをインストールするの？」というところからつまずきがちです。
そんな時は、以下の書籍をおすすめします。
改訂新版 OpenFOAMの歩き方 (技術の泉シリーズ（NextPublishing）)
インストール方法とチュートリアルで流体解析を体験・・・ちょっと高度な解析まで解説があります。著者曰くOpenFOAMのバージョンの追跡を行いながら、書籍をアップデートするようなので安心ですね。
Follow @t_kun_kamakiri
