第38回Modelicaライブラリ勉強会「FMU4FOAMによる竹とんぼシミュレーション（その２）」

2023年1月29日 2023年12月18日

openscc

　OpenFOAMとOpenModeliacaのFMUを連成するFMU4FOAMライブラリを用いた竹とんぼシミュレーションに取り組んでいます。内容の理解が進みOpenFOAMについてもOverSetMeshを使うともっと簡単にモデル作成ができることが分かりました。
　そして、前はできなかった竹とんぼが上方に飛行する３Dアニメーションもできることがわかりました。FMU4FOAMについては以前の記事などをご参照ください。
　第38回Modelicaライブラリ勉強会に参加して新たなモデルでのシミュレーションの方法について情報提供しました。

　以前のモデルについて説明します。OpenFOAMは竹とんぼの空力シミュレーションを行います。ソルバーはpimpleFoamを使いました。OpenModelicaからのz軸回転角度angle_outと回転中心CofR_out(x,y,z)を使って空気中の赤の領域のメッシュを回転し、竹とんぼにかかる流体力とトルクを計算します。物体は回転運動のみで並進運動はしません。メッシュを移動させるz軸の回転角度を時間的に変化させてメッシュを回転させるライブラリはなかったので新たに開発しました。
　OpenModelicaでは、剛体運動用のbodyとz軸回転運動用のinertiaで構成されています。inertiaではトルクt_inを入力して回転運動の計算を行いz軸回転角度angle_outと角速度omega_outを出力します。bodyでは力f_inとt_inを入力して剛体運動の計算を行い、重心pos_outを出力します。重心と角速度はOpenFOAMでは使いません。
　本来、bodyのみで6自由度の並進と回転の剛体運動シミュレーションが可能ですが、z軸のみの角度と回転速度の出力方法がよくわからなかったため少し複雑なモデルになっています。

　次に、今回の竹とんぼ飛行シミュレーションモデルについて説明します。
　OpenFOAMではoverSetメッシュを使いました。backgroudメッシュとoverSetメッシュと２つのメッシュを用意して、overSetメッシュの境界でbackgroundメッシュの値を補間する方法です。メッシュ変形がなく、計算の不安定化や精度劣化が少ないため移動物体などのシミュレーションによく利用されています。ソルバーはoverPimpleDyMFoamを使いました。
　OpenModelicaから６自由度の重心pos_in、回転中心CofR_in、回転角度RPY_inが入力されます。回転中心CofR_inは重心を使用します。これらからoverSetメッシュを６自由度移動します。前回は１軸の回転運動だけの移動でしたが、３自由度の並進運動、３自由度の回転運動が計算できます。物体移動に伴って流体運動が生じる力f_outとトルクm_outを計算します。
　OpenModelicaは基本的にbodyで構成されています。OpenFOAMから流体運動による力f_inとトルクt_inが入力されて、剛体運動シミュレーションを行い重心pos_outと回転角度angle_outを計算します。前回に比べてモデルが簡単ですっきりしています。

　計算モデルについて少し詳しく説明します。backgroundメッシュの領域は600mm角の立方体です。overSetメッシュは直径160ｍｍ高さ190mmの円柱でその中に竹とんぼの形状が入っています。overSetメッシュの作成方法についてはこちらをご覧ください。
　層流モデルで乱流は考慮していません。動粘性係数は空気のものを使用しています。overSetメッシュを利用する場合は、dynamicFvMeshライブラリとしてdynamicOversetFvMeshを使います。overSetメッシュの移動方法を指定するsolidBodyMotionFunctionとして、並進運動用のcoupledTranslationMotionと回転運動用のcoupledRotationMotionを使います。FMU4FOAMの標準ライブラリです。
　初期条件としてz軸回転速度を500rad/sとしました。
　計算ステップdeltaTimeは1e-4秒と固定です。OpenModelicaも同一にしています。OpenFOAMでは最大クーラン数maxCoなど使って計算状況に応じて計算ステップを可変にできますが、今回の計算では計算が不安定になるので固定間隔としています。

　計算メッシュを示します。69065セルです。

　dynamicFvMeshDictです。overSetメッシュだけをメッシュ移動の対象とするために、cellSetでmovingZoneを指定しています。並進運動ではcoupledTranslationMotionを使い、移動位置名posNameをpos_inとしています。回転運動はcoupledRotaionMotionを使い、回転角度名rollPitchYawNameをRPY_in、回転中心名CofRNameをCofR_inとしています。

　system/externalCouplingDictです。OpenFOAMから出力すｒ力とトルクに関する設定をします。typeはextForcesを入力します。力とトルクを計算する物体のpatch名patchesを"taketonbo"と指定しています。流体密度rhoはrhoInf=1.2kg/m3としています。力の名前forceNameをf_out、トルクの名前momentNameをm_outとしています。

　controlDictです。FMU4FOAMライブラリはexternalCommでlibsでインクルードする必要があります。ただ、並列計算する場合に領域分割のためdecomposeParコマンドを使いますが、エラーとなります。コメントアウトすると領域分割できます。その後、もとに戻して計算を実行することができます。
　FMU4FOAMを使うには、functionsにFMUSimulatorを指定する必要があります。

　OpenModelicaのモデルを詳しく見てみます。
MultiBodyではworldが必須です。座標と重力方向を設定します。基本的に３次元６自由度の剛体運動の標準モデルbodyを使います。質量、慣性モーメントなどを設定します。
　OpenFOAMからの入力として３成分トルクt_inと３成分力f_inを追加します。3成分の入力とするには、Attributesを開いて、Dimensionsに[3]と入力する必要があります。
　力とトルクの両方をbodyに作用するモデルとしてforceAndTorqueを使用します。直線矢印に力入力を接続し、曲がった矢印にトルク入力を接続します。
　bodyの計算結果をOpenFOAMへ外部出力するためにabsoluteSensorを使います。位置、速度、加速度、角度、角速度、角加速度の６つの出力が可能であります。矢印の位置で出力が決まっています。左から一番目の３成分位置pos_outと４番目の３成分角度angle_outを使います。

worldの設定です。label2に"z"を入力し、nを{0, 0, -1}とします。

　bodyの設定画面を説明します。r_CMに{0, 0, -0.0253}と重心位置を設定します。mは質量で0.00494kgを入力します。慣性モーメントは６成分ありますが、ここではz軸まわりのI33成分に0.0451kgm2を入力します。

absoluteSensorの設定です。出力したいものにチェックします。今回はget_rとget_anglesです。

　計算結果を示します。60msまでの竹とんぼが回転しながら上方に移動している様子がシミュレーションできました。計算時間は822秒です。重心位置position、物体角度angle、力force、トルクtorqueの履歴をグラフに示します。時間の経過にともない、竹とんぼがz軸回りに回転しながらz方向に移動し、z軸回りのトルクが作用するので、回転速度が減速して揚力が低下していくことがわかります。
　OpenModelicaでは角度が一回転で-3.14~3.14の範囲で変わりますが、OpenFOAMでは連続して解釈して計算が実行できるようです。但し、一回転ごとに力やトルクで不連続な部分が見られます。