私と円柱と抵抗と - 人生の人生

この記事は東京大学航空宇宙工学科/専攻 Advent Calendar 2018の21日目の記事です．

0. はじめに
1. 解析対象
2. 計算条件
2.1. 問題設定等
2.1. CFDについて
3. 計算結果
4. まとめ
参考文献

0. はじめに

どうもこんばんは，人生です．主な生息場所はTwitterです．Twitter社さん，もう凍結しないでください．
普段は所属研究室内製のCFDソルバで遊んでいたりします．
本稿でもCFDソルバで遊んでみた結果を書きます．

1. 解析対象

本稿では，空気抵抗を受ける2次元円柱の弾道計算を行います．
本当は球を用いた3次元解析をして現実世界での結果と比較～～～！なんてことをしたいのですが，計算コストがエグいので2次元円柱 (つまり，円) を用いた2次元解析を行います．
円柱や球に働く空気抵抗は実験値や近似値，解析値が出回っていますが，本稿ではこの空気抵抗値を所属研究室内製のCFDソルバで推算します．
特に，円柱の運動とともに対気速度やレイノルズ数が変化すると考え，その刻々に対して空気抵抗値を更新しながら運動を発展させることを狙います．
結局のところただのルンゲクッタなので読む価値はないです．
ていうか，空気力推算でミスってる気しかしないです．
レイノルズ数ってなんだっけ...．無次元化ってなんだっけ...．

2. 計算条件

2.1. 問題設定等

円のサイズは野球の硬式球を想定して直径74 mm，質量145 gとします．
これに初速30 m/sを持たせて，30 degの投射角で0 ℃，1 atmの静止気体中に投げます．

また，下図のように座標系と力を定義します．

f:id:sktakefumi:20181213232104p:plain:w300 — 座標系，力の定義

$L$ , $D$ はいわゆる揚力，抗力です．
本稿では静止気体中の運動を考えているので図中の向きとなります．
運動方程式は以下のようになります．
$\left\{ \begin{align} \dot{x}(t) &= u(t) \\ \dot{u}(t) &= F_x(t)/m \\ \dot{y}(t) &= v(t) \\ \dot{v}(t) &= F_y(t)/m - g \\ \end{align} \right.$
さらに， $L$ , $D$ は揚力係数 $C_l$ ，抗力係数 $C_d$ を用いて次のように書かれます．
$\left\{ \begin{align} L &= \frac{1}{2} \rho |\vec{V}|^2 S C_l \\ D &= \frac{1}{2} \rho |\vec{V}|^2 S C_d \\ F_x &= - L \sin \gamma - D \cos \gamma \\ F_y &= L \cos \gamma - D \sin \gamma \\ \end{align} \right.$
力の定義の仕方に頭の悪さがにじみ出ていますが，ソースコードではちゃんと
$\begin{align} \dot{v}(t) &= F_y(t)/m \\ F_y &= L \cos \gamma - D \sin \gamma - mg \\ \end{align}$
としているので許してください (?)

再三になりますが，本稿では $C_l$ ， $C_d$ をCFDソルバを用いて推算します．
なお，この流れでは $C_l=0$ となるはずなので，0以外の値は計算上のミスとして0で上書きします．
じゃあなんで定義したんだ，人生くん...．
また，CFD解析は各時間ステップで初期状態から計算するため，前の時間ステップで形成された流れ場は次の時間ステップには影響を及ぼさないものとします．
これは，CFDソルバが動体解析に対応していないためです．
先日出席した学会ではいろんな人が平然と非定常構造連成解析をしててビビりました．

後は普通に4次のルンゲクッタにぶち込みます．
この際， $t+dt/2$ の時点でも $C_l$ ， $C_d$ を推算しようとすると計算量がエグいことになるので， $t$ – $t+dt$ 間の $C_d$ は $t$ における値で一定とします．

2.1. CFDについて

本稿では，支配方程式をRANS方程式として，2次元定常計算を行います．
計算スキーム等は省きますので，詳しいことやそのほかソルバについては参考文献[1]をご覧ください．

計算格子に関して，本稿では申し訳程度の研究 () 要素として，解適合格子細分化法 (Solution-Adaptive Mesh Refinement) を利用します．
最近はもっぱらこれをいろんな解析対象に適用して計算結果の向上を目指すみたいな，小学生の自由研究にも満たないことをやっています．
頭が欲しいんだ．

本解析では例えば下図のような格子が生成されます．

f:id:sktakefumi:20181220215254p:plain — AMRによる計算格子の発展

これにより，円柱からの後流が，初期格子では素早く拡散してしまうのに対し，AMR格子では風下方向に長く維持されるようになります．

f:id:sktakefumi:20181220215726p:plain — マッハ数 $M$ 分布，気流は左から右，左: 初期格子 (4,356セル)，右: 3サイクル適合格子 (8,376セル)

f:id:sktakefumi:20181220215733p:plain — マッハ数 $M$ 分布，気流は左から右，左: 初期格子 (4,356セル)，右: 3サイクル適合格子 (8,376セル)

この解像度の向上は $C_d$ に若干の感度を有し，初期格子の234 cnt (カウント) からAMR格子の226 cntまで，3.4 %ほどの変化量をもたらします．
あんま意味ねえな...．

f:id:sktakefumi:20181220220739p:plain — 抵抗係数 $C_d$ の収束履歴，擾乱位置でAMRを適用している

ただ正直なところ，抵抗係数については最小格子幅 $\Delta x$ が支配するところも大きく，セル数の増分を壁面と空間のどちらに割くかというのは微妙なところです．
現時点では，摩擦抵抗 $C_{d,f}$ を突き詰めたければ壁面の，圧力抵抗 $C_{d,p}$ を突き詰めたければ空間の解像度を上げる必要があるという見解です．

3. 計算結果

こんな感じになりました．

f:id:sktakefumi:20181220222631p:plain — 軌道の計算結果

ただただ,「それっぽい」絵が描けたというだけでした．
個人的にはそんなに前に落ちるのかというのが驚きですが，計算をミスってる感がどうもぬぐえないのでそういうことかもしれません．
本当はエネルギーの収支とかを考えたら面白いのかもしれませんがなんかあまりにもしょうもなさすぎてやる気を失ってしまったのでこれで終わります．

4. まとめ

まじでしょうもねえ．
助けてください．

参考文献

[1] Yoshiharu Tamaki and Taro Imamura. "Turbulent Flow Simulations of the Common Research Model Using Immersed Boundary Method", AIAA Journal, Vol. 56, No. 6 (2018), pp. 2271–2282. https://doi.org/10.2514/1.J056654