2019年4月10日水曜日

pimpleDyMFoamを用いた流体と剛体運動の連成解析


はじめに

運動している剛体周りの流体挙動とその連成を数値解析的手法で調べる。
OpenFoamのソルバーにpimpleDyMFoamがある。移動格子での流体解析が実装されており、これによって剛体周りの流体とその運動を記述できるかもとのこと。
ここでは１軸自由度の流体運動連成解析を考える、剛体運動は以下の運動方程式で記述する。
\begin{align} I\frac{\partial^2 \theta}{\partial t^2} = M \end{align}
流体によるモーメントから上記方程式に従って迎角が得られる、それを移動格子を通じて流体解析にフィードバックする。格子移動時の物理量内挿はALEでやってる？

以下、いろいろ適当な点/間違いあるかと思います。

計算対象

なんでも良いのだけどGrabCADから適当なumbrellaのSTLを見つけてくる（傘の運動に興味があった）
umbrella.stl

適当にリスケールしています、だいたい代表長0.8mくらい。

計算条件

ヴァージョン：OpenFoam v5.0
サーバーOS：Ubuntu Server 16.04.6 LTS
並列数：72

ベースにしたチュートリアルケース：/opt/openfoam/openfoam5/openfoam5/tutorials/incompressible/pimpleDyMFoam/wingMotion
（インストール場所/opt/openfoam/openfoam5です）

まず次のようなスクリプトを準備
./qfoam_header.sh 
#!/bin/bash

ncpu=72
dir_pre_mesh=mesh
#dir_main_simple=simpleFoam
dir_main_dym=pimpleDyMFoam

source /opt/openfoam/openfoam5/openfoam5/etc/bashrc

# Source tutorial run functions
. $WM_PROJECT_DIR/bin/tools/RunFunctions

計算格子の作成

計算格子はsnappyHexMeshで作成する。
流体解析的な精度は今回は気にしないので格子は真面目に作らない。200万cells数くらいと思う。

以下を編集
{dir_pre_mesh}/system/controlDict
{dir_pre_mesh}/system/blockMeshDict
{dir_pre_mesh}/system/snappyHexMeshDict
{dir_pre_mesh}/system/surfaceFeatureExtractDict
{dir_pre_mesh}/system/decomposeParDict

#!/bin/bash

. qfoam_header.sh

# Make mesh
cd ${dir_pre_mesh}

cp -r 0.orig 0

surfaceFeatureExtract > log_surfaceFeatureExtract 2>&1
blockMesh > log_blockMesh 2>&1
decomposePar > log_decomposePar_1 2>&1
mpirun -np $ncpu snappyHexMesh -overwrite -parallel > log_snappyHexMesh 2>&1
#snappyHexMesh -overwrite > log_snappyHexMesh 2>&1

# For parallel process
reconstructParMesh -constant > log_reconstructParMesh 2>&1
rm -r processor*
decomposePar > log_decomposePar_2 2>&1

計算格子を作成した。

なぜか傘の柄の一部が消えてしまった。。

格子ファイルのコピー

メイン計算は別ディレクトリで行う。
#!/bin/bash

. qfoam_header.sh

# Copy the mesh
cd ${dir_main_dym}
cp -r ../${dir_pre_mesh}/processor* ./
cp -r ../${dir_pre_mesh}/0 ./
cp 0.orig/pointDisplacement ./0/
cp -r ../${dir_pre_mesh}/constant/polyMesh ./constant/
#mapFields ../${dir_pre_mesh} -parallelTarget -consistent -sourceTime latestTime -parallelSource > log_mapFields 2>&1
#
count=0
while [ $count -lt $ncpu ] ;
do
  mv processor$count/0/pointDisplacement.unmapped* processor$count/0/pointDisplacement
  count=`expr $count + 1`
done

物理モデル・計算モデル

分子粘性係数：1.e-5 Pa.s
乱流モデル：なし（層流計算）
一様流速度：10 m/s
一様流密度：1 kg/m3 (メイン計算では使わないけど空気力算出に必要)
Re数はだいたい1e6

以下を編集
{dir_main_dym}/constant/transportProperties
{dir_main_dym}/constant/turbulenceProperties
{dir_main_dym}/constant/dynamicMeshDict
慣性モーメントや重心、質量はすごく適当な値になっている。

流体＋移動格子計算

次を編集
{dir_main_dym}/system/controlDict
{dir_main_dym}/system/decomposeParDictの"numberOfSubdomains 72;"も編集する必要があるはず。

pimpleDyMFoaMソルバーを実行する
#!/bin/bash

. qfoam_header.sh

# Solve flowfield
cd ${dir_main_dym}
mpirun -np $ncpu pimpleDyMFoam -parallel > log_pimpleDyMFoam 2>&1
しばらく待ちます・・5秒間の計算に実時間50時間程度かかった。

可視化

計算結果からVTKファイルを作成する
#!/bin/bash

. qfoam_header.sh

cd ${dir_main_dym}

#reconstructPar -latestTime > log_reconstructPar_post 2>&1
reconstructPar > log_reconstructPar_post 2>&1
foamToVTK > log_foamToVTK 2>&1

ParaView 5.6.0で可視化
・t=3.0sのときの速度ベクトル分布


空気力と迎角の時間履歴

ピッチングモーメント・迎角いずれも収束に向かっており、動的不安定性は出なかった。

おわりに

計算モデルとしてそれなりにロバストな印象はあるが、もう少し高速化できればいいかもしれない。
なお、いろいろ調べていくとやっぱり乱流モデルの影響は無視できないっぽい。

参考文献等

参考にさせていただきました：

投稿者 時刻:

0 件のコメント:

コメントを投稿

登録: コメントの投稿 (Atom)