「テーマ」

・OpenFoam2D円柱周りのCd解析結果を実験値と比較する。

「方針」

・幾つかのレイノルズ数において2D円柱周りの乱流流体解析を実施し、風洞試験値と比較してみます。

　尚Openfoamを用いた2D翼の抵抗係数確認にて幾つかの乱流モデルの比較を行いました。

結果、pisoFoam^*1+RAS^*3(k-epsilon)よりもsimpleFoam^*2+RAS(k-OmegaSST)の方が風洞試験値と一致することを確認しています。それを受けて本稿でもsimpleFoam+RAS(k-OmegaSST)にて計算してみます(ソルバー詳細：参考文献[1]2.1項)。

　※1 pisoFoam：非圧縮・非定常乱流解析ソルバー(PISO法)、時間経過を観察できる。

　※2 simpleFoam：非圧縮・定常乱流解析ソルバー(SIMPLE法)、定常状態(収束状態)を求める。

　※3 RAS：Reynolds Averaged Simulationの略。速度を平均値とそこからの偏分という形にしてナビエストークスの式を展開する。この際用いられる乱流モデルには幾つかの種類があり、代表的なものとして以下を選択した。

　k-epsilonモデル：流体解析で最初に出てくる標準モデル

　k-OmegaSSTモデル：境界層の振舞いを改善したモデル

---

１．１　　解析モデル＆OpenFoam設定

　円柱回りの解析を行うためのモデルを設定します。

・空間メッシュ

・円柱直径1[m]。但し計算は半円(真円では上下流が干渉を起こして収束しないと推定)。

・メッシュ：16(長さ) x 3(高さ) x0.05(幅)[m]を160x30x1の合計4800マスに分割。

　円柱周辺を4、8倍と段々に細分化、円柱表面には約1[mm]厚を30層設定。

・周囲速度0.01、0.1、1、10、30[m/s]の5ケース、上下面はslip条件、円柱表面をnoSlip条件。

・圧力：出口をゲージ圧0で固定

(この方が安定します。但し円柱後流領域を十分取る必要があります)。

・空気動粘性：1.5e-5 [m²/s]

・抵抗係数Cdの計算は追加機能関数forceCoeffsを使用。

・解析ソルバー：simpleFoam(非圧縮・定常乱流解析ソルバー)

・乱流モデル：k-OmegaSST

・k及びomegaの初期値は適当(I=1%で参考文献[1]4.5.5項に従う。)

・nutの初期値は0。

上記以外のパラメータは以下のOpenFoamチュートリアルの物をそのまま使用。

  \tutorials\incompressible\simpleFoam\motorBike

　また比較対象となる風洞試験値として1921年Wieselsbergerによって示された円柱のレイノルズ数別Cd値の実験値を選びました(原図からの読み取りのため概略値)。

---

１．２　OpenFoam計算

　ケース別に計算結果を以下に示します。

〇ケース１：周囲速度0.01[m/s](円柱頂上Rex：5e+2)

結果、Cd=1.05に収束しました(半円を円とするため0.525を2倍)。これはWieselsberger試験値1.2(Re<3e+5)に大分近い値を示しています。また剥離点は円柱頂上より前方側(よどみ点より75度位の所)に確認出来ます。試験値(参考文献[2]9.2項)では80度位とのことでおよそ一致しています。剥離が早い分円柱表面では主流の圧力回復が進まず、抵抗が大きくなる様子が伺えます。

　ちなみにこの状態でy+は円柱前方で0.02、剥離点近傍で0.12と十分なメッシュ細かさとなっています(境界層内部の速度分布を観察するためy+は壁関数に頼らない1以下が望ましいと考えました)。また円柱頂上での推定境界層厚さ^*4は83[mm]に対しメッシュ厚さが1[mm]なので十分細かいメッシュ層と考えます。

　*4：推定境界層厚さ　乱流境界層と想定し、以下の経験式より算出(参考文献[2]4.8.1項)。

　　δ=0.37*x/Rex^0.2

　　Rex=U*x/ν

　　U：周囲速度0.01[m/s]

　　x：代表長さ　0.5*2π/4[m]

　　ν：空気動粘性1.5e-5[m²/s]

〇ケース２：周囲速度0.1[m/s](円柱頂上Rex：5e+3)

結果、ケース１とほぼ同じCd=1.02に収束しました。また剥離点はケース１とほぼ同じ辺りに確認出来き、実験値(参考文献[2]9.2項)80度位とおよそ一致しています。

　ちなみにこの状態でy+は円柱前方で0.11、剥離点近傍で0.67と十分なメッシュ細かさとなっています。また円柱頂上での推定境界層厚さは52[mm]に対しメッシュ厚さが1[mm]なので十分細かいメッシュ層と考えます。

〇ケース３：周囲速度1[m/s](円柱頂上Rex：5e+4)

　解析の結果、Cdが振動し続けます。これは解析対象が剥離というダイナミックな運動のためと考えます。収束が望めなさそうなので30000stepのうち、1000step毎の値の平均をとりました。結果、Cd=0.636となしました。これはWieselsberger試験値の1.2(Re<3e+5)と0.35(Re>3e+5)のちょうど中間を示しています。図より剥離点は円柱頂上辺りに確認出来ます。剥離が遅れた分、円柱表面では主流の圧力回復がケース１より進んだため、抵抗が小さくなっている様子が伺えます。

　ちなみにこの状態でy+は円柱前方で0.58、剥離点近傍で4.4とそれなりのメッシュ細かさとなっています。また円柱頂上での推定境界層厚さは33[mm]に対しメッシュ厚さが1[mm]なので十分細かいメッシュ層と考えます。

〇ケース４：周囲速度10[m/s](円柱頂上Rex：5e+5)

　こちらもCdが振動を続けたため30000stepのうち、1000step毎の値の平均をとりました。結果、Cd=0.372となしました。これはWieselsberger試験値0.35(Re>3e+5)と良く一致しています。また剥離点は円柱頂上より後方側(よどみ点より110度位の所)に確認出来ます。実験(参考文献[2]9.2項)では120度位とのことで良く一致しています。注目すべきは110度位で一度剥がれた気流が120度位のところに再付着しようとしている流れの様子が見て取れます(完全に再付着はしていませんが)。

　ちなみにこの状態でy+は円柱前方で2.8、剥離点近傍で32と若干メッシュ細かさが不十分に見えます。一方円柱頂上での推定境界層厚さは21[mm]に対しメッシュ厚さが1[mm]なので厚さ方向は十分細かいメッシュ層と考えます。

〇ケース５：周囲速度30[m/s](円柱頂上Rex：1.6e+6)

　Cdは振動を続けたため30000stepのうち、1000step毎の値の平均をとりました。結果、Cd=0.396となしました。これはWieselsberger試験値0.35(Re>3e+5)と良く一致しています。また剥離点は円柱頂上より後方側(よどみ点より110度位の所)に確認出来、実験値120度位と良く一致しています。注目すべきはこちらも10[m/s]と同様に110度位で一度剥がれた気流が120度位のところに再付着しようとしている流れの様子が見て取れます。

　ちなみにこの状態でy+は円柱前方で5.2、剥離点近傍で92と相当メッシュ細かさが不十分となっています。一方円柱頂上での推定境界層厚さは17[mm]に対しメッシュ厚さが1[mm]なので厚さ方向は十分細かいメッシュ層と考えます

---

２．まとめ

　上述計算結果をWieselsberger試験値グラフに重ねてみます。

　　(プロットにRex(局地レイノルズ数)を用いたため試験値Re(円直径)と若干ずれます。)

　Rex=5e+4(1[m/s])にて試験値からずれていますが、それ以外は試験値とよい一致を示しています。

　本計算は乱流境界層剥離点の移動を伴うダイナミックな事象が対象であり、試験値のような急変を模擬することは難しいと思えるため、上記結果は十分な模擬度と感じます。

　この計算結果をまとめてみます。

・臨界レイノルズ数(3e+5)以下では剥離点が円柱頂点よりも前方にあり、主流の圧力回復が阻害される。

　結果抵抗係数Cdは大きくなる。この様子がOpenfoamの流線からも確認出来た。

・臨界レイノルズ数(3e+5)以上では剥離点が円柱頂点よりも後方にあり、主流の圧力回復が進む。　　

　結果抵抗係数Cdが小さくなる。Openfoamの流線から剥離後の渦による再付着の様子も見て取れる。

注意：今回のk-OmegaSSTモデルでは初期境界条件として以下の壁関数を用いています。それにも関わらすメッシュをy+の適正値よりも遥かに小さくとっているため、そこの挙動に異常がないかは常日頃注意する必要があると考えます(今回は異常な挙動をしていません)。

　　nut：nutkWallFunction

　　k：kqRWallFunction

　　Omega：omegaWallFunction

　尚今回の乱流モデルkOmegaSSTは全領域を乱流境界層として計算を行います。しかし円柱前方のよどみ点辺りは層流境界層と予想されます。そこで「補足１」に乱流モデルを層流モデルにした場合の計算結果を示します。結果として層流モデルは試験値の違いが大きく、ある程度のレイノルズ数(1e+2)からは乱流モデル(k-OmegaSST)が妥当と考えます。

　また今回の本文の計算では予想される境界層よりも十分細いメッシュを用意しました。しかしこの精細度のメッシュを3Dで用意することはPC(パソコン)解析では厳しいものがあります。そこで敢えて本文よりもメッシュを粗くした計算を「補足2」にて実施してみます。結果、メッシュを境界層厚さ程度まで粗くしても計算結果に大きな違いは出なさそうです。

　更に「補足3」にpisoFaomを用いた場合の注意点を参考として示しておきます。

---

「補足１」

　本文ではk-OmegaSST乱流モデルを使用したが、円柱よどみ点付近は層流境界層が予想されます。そこで本文計算と同じ条件で層流モデルの場合を

計算してみます。

　変更点は以下のみとなります。

・turbulencePropertiesファイル内構文を変更

　(新)simulationType 　　laminar;

　(旧)simulationType      RAS; RAS{RASModel  kOmegaSST;  turbulence on;}　結果は以下となりました(速度0.1[m/s]を実施)。

Cd=0.602で乱流モデル1.02よりも大分小さくなっており、試験値からも離れています。

この結果をもう少し詳しく見てみます。

　抵抗係数Cdを前後圧力差による形状抵抗Cd_pressと摩擦による摩擦抵抗Cd_viscosに分けて考えます。paraviewより上記計算結果を分解すると以下になりました。

　乱流モデル　1.02(Cd_RAS)     = 0.99(Cd_press)+0.03(Cd_viscos)

　層流モデル　0.602(Cd_laminar) = 0.58(Cd_press)+0.022(Cd_viscos)

　この結果より摩擦抵抗はそれ程違いがないのに対し(それでも乱流の方が大きいです)、形状抵抗に大きな違いが出ていることが分かります。

　これは主流の様子が大きく違うことを意味していますが、全体流れの図を見ると確かに流れの様子が大きく違います。興味深い点は層流モデルの方が円柱の後方で大きく主流が跳ね上がっている点です。

　いずれにしても本比較より円柱よどみ点辺りは層流境界層だと考えられますが、それ以降はかなり早い時点で乱流境界層に切り替わっていると想像されます。

---

「補足２」

　本文計算のメッシュを敢えて粗くした計算を行ってみます。

　まずはメッシュ厚さを想定される境界層と同レベルとした計算を実施します。

　・円柱直径1[m]。但し計算は半円。

　・メッシュ：16(長さ) x 3(高さ) x0.05(幅)[m]を160x30x1の合計4800マスに分割。

　　円柱周辺を4倍に細分化、円柱表面には約20[mm]厚を3層設定。

　他の条件は1.1と全く同一とします。

　計算結果をWieselsberger試験値グラフに重ねてみます(今回の計算は粗メッシュとして表示)。

　精細メッシュと比べて若干試験値から離れていますが、それなりの結果を示しています。

　続いてメッシュ厚さを想定される境界層の2.5倍のレベルとした計算を実施します。

　・円柱直径1[m]。但し計算は半円。

　・メッシュ：16(長さ) x 3(高さ) x0.05(幅)[m]を160x30x1の合計4800マスに分割。

　　円柱周辺を2倍に細分化、円柱表面には約50[mm]厚を2層設定。

　他の条件は1.1と全く同一とします。

　計算結果をWieselsberger試験値グラフに重ねてみます(今回の計算は超粗メッシュとして表示)。

　ここまでメッシュを厚くするとかなり試験値から離れてしまいます(超臨界前後の平均に近いようです)。

　この変化をもう少し詳しく見てみます。メッシュを粗くすると大きく影響の出る一つとしてnut(乱流動粘性)があります。

OpenfoamのHPよりnutの壁関数は下式で定義されています。

　　ν_t：乱流動粘性

　　ν：空気動粘性

　　ｙ+：無次元壁面距離

　　E：対数則の積分定数に相当する経験定数。(通常9.8:参考文献[1]5.7.9項)

　　κ：カルマン定数。(通常0.4～0.45:参考文献[1]5.7.9項)

　そこで周囲速度0.1[m/s]でのnutの様子を精細メッシュと超粗メッシュで比較してみます。

　精細メッシュではnutが成長しだすのが、円柱頂上の後方なのに対し、超粗メッシュでは円柱頂上前方で成長しだしています(流れに関係なく壁関数に従ってしまう)。nutの増加は動粘性の増加を意味しており、動粘性の増加によって境界層上下のせん断力が増加します。これにより境界層上下で運動量の伝搬が起こり、結果剥離が遅れると考えます。剥離が遅いと主流の圧力回復が促進され、抵抗の減少につながります(この説明、風が吹けば桶屋が儲かるに近い説明になっており、あくまで推測です)。

　以上より可能ならばメッシュ・サイズアップはせいぜい想定される境界層厚さまでとなりそうです。

---

「補足３」

　ここではpisoFoam+k-epsilonの組み合わせにおける失敗例を示しておきます。

〇解析モデル＆OpenFoam設定

　下記のような2D円柱を設定します。

・円柱直径1m

・メッシュ：16(長さ) x 6(高さ) x0.04(幅)[m]を160x60x1の合計9600マスに分割。

　円柱周辺の1.5(長さ)ｘ1.5[ｍ](高さ)の区間を2倍に細分化。

・周囲速度0.01、0.1、1、5、10[m/s]の５ケース、上下面はslip条件、円柱表面をnoSlip条件。

・圧力：出口をゲージ圧0で固定

・空気動粘性：2.0e-5 [m²/s]

・解析ソルバーはpisoFoam(非圧縮・非定常乱流解析ソルバー)

・乱流モデルはRAS(レイノルズ平均に基づく乱流モデル)でkEpsilonモデル(標準ｋ-εモデル)

・epsilon、k及びomegaの初期値は適当(I=1%で参考文献[4]4.5.5項に従う。)

・nutの初期値は0とする。

上記以外のパラメータは以下のOpenFoamチュートリアルの物をそのまま使用。

  tutorials\incompressible\pisoFoam\RAS\cavity

　結果出来上がったメッシュは円柱表面で厚さ50[mm]。10[m/s]で想定される乱流境界層厚さは1.2項ケース４より21[mm]で境界層内部の流体挙動を観察するにはかなり足らない状態となります。これは「補足２」での後半の計算とほぼ条件が一致しており、simpleFoam+k-OmegaSSTでは抵抗係数Cd=0.6位に落ち着いていました。

〇計算結果(クーラン数0.1)

　クーラン数Coが0.1となるように積分幅を設定して計算した結果、Cd値は以下となりました。

　　Cd=1.21(周囲速度0.01[m/s])@Rex=5e+2

　　Cd=1.05(周囲速度0.1[m/s])@Rex=5e+3

　　Cd=1.02(周囲速度1[m/s])@Rex=5e+4

　　Cd=1.01(周囲速度5[m/s])@Rex=2.5e+5

　　Cd=0.98(周囲速度10[m/s])@Rex=5e+5

　試験値に重ね合わせるとそれなりに一致しているようです。但し超臨界(Re＞3e+5)の挙動は再現できていません。

〇計算結果(クーラン数0.01)

　続いてクーラン数Coが0.01となるように積分幅を設定して計算した結果、Cd値は以下となりました。

　　Cd=1.17(周囲速度0.01[m/s])@Rex=5e+2

　　Cd=0.34(周囲速度0.1[m/s])@Rex=5e+3

　　Cd=0.27(周囲速度1[m/s])@Rex=5e+4

　　Cd=0.25(周囲速度5[m/s])@Rex=2.5e+5

　　Cd=0.26(周囲速度10[m/s])@Rex=5e+5

　試験値に合わせると今度は超臨界(Re＞3e+5)の挙動を再現しているように見えます。

　しかしこのようにクーラン数で計算結果が変わってしまってはどちらを信用すべきか判断がつきません。

　この失敗より以下の３点を注意する必要を反省しました。

・境界層を議論する場合はメッシュ粗さが境界層厚さよりも小さいことが望まれる。

・pisoFoamは動的挙動を扱うためクーラン数による計算結果の違いが発生しうる(pisoFaomは反復計算が少ない)。

・乱流モデルの設定に注意。乱流境界層を扱うのであればk-epsilonでは役不足。

---