OpenFOAMによる球体の抗力係数調査
-
製品名
- OpenFOAM
-
カテゴリー
- 熱流体解析事例
-
産業分野別
1．はじめに
自動車やバイク等の設計において、抗力係数は製品性能に直結する重要な指標となります。また、ビルなどの建造物設計においても、強度設計のみならず風環境評価の観点から欠かせない要素です。
抗力係数は風洞実験で計測することも可能ですが、流体解析（CFD）を活用することで、形状変更への迅速な対応や実機スケールでの性能評価が可能となり、開発コストや期間を大幅に削減できます。そのため、抗力係数を精度よく予測（計算）できるかは、流体解析の妥当性や信頼性を評価する上で極めて重要です。
球体の抗力係数は、レイノルズ数1.0×103～2.5×105の範囲では0.4～0.44の値を、また臨界レイノルズ数は2.5×105～3.5×105の範囲に位置することが実験により知られています。
本解析事例では、Foundation版 OpenFOAM ver8の非定常-非圧縮の単相流ソルバーpimpleFoamを利用して、乱流場における球体の抗力係数を算出し、公開されている実験データと比較しました。
2．モデル説明
図1に、計算モデルの概要図を示します。 矩形の計算領域に球体を配置し、一定の速度で空気を流入させます。 球体の直径は0.3 mで、乱流モデルはk-ω SSTとしました。 また、流入境界に設定した流速は、 Allen域からNewton域の移行区間から極超臨界領域までの0.1 、1、10、50、100、200 km/hの6ケースとしました。
表1に計算ケースとレイノルズ数を示します。
-
図1．モデル概要
表1．計算ケース
補足①：乱流モデルk-ω SSTについて
計算に使用した乱流モデルk-ω SSTについて説明します。
k-ω SST はレイノルズ平均ナビエストークス式(RANS)の2方程式モデルで、標準 k-ε と k-ω のハイブリッドの乱流モデルです。
RANSのベースとなる標準k-εモデルは、乱流エネルギーkと散逸率εの2つの輸送方程式をそれぞれ解くことで、乱流場を解く計算手法です。流体解析で広く利用されてきたモデルですが、問題点として、壁面近傍の計算精度が低いことが知られています。y+も20～200程度の範囲に収まるように、壁面近傍の計算点を配置する必要があります。
一方、k-ωは、壁面近傍の予測精度に優れた乱流モデルで、乱流エネルギーkと乱れの周期ω(＝ε/k)の、2つの輸送方程式を解く手法です。壁面近傍の精度がk-εより高く、y+も1（＝粘性低層）近くに設定することが可能です。なお、流入条件としてωを0とすると、渦粘性μt(=ρk/ω）が0割で不定となるため、有限で微小な値を設定する必要があります。
k-ω SSTは、標準k-εとk-ωの欠点を補うハイブリッドの乱流モデルで、壁面近傍ではk-ωを使って計算し、壁面から十分離れた完全乱流域ではk-εに切り替える計算手法です。 3次元-非定常の空力計算で乱流モデルを使用する際は、選択肢としてより高精度なLESもありますが、工業用途では計算コストと精度のバランスに優れることから、
k-ω SSTが広く利用されています。
3．計算条件
図2に、境界条件を示します。計算領域の大きさは、球によって生じる乱れが境界の影響を受けないように、境界の位置を球から十分離れた位置に配置しました（流路断面積に対する球の断面積の比率は4％以下）。また、全てのケースで、流入境界における乱流強度を、0％相当の十分小さい値としました。
-
図2．境界条件
図3に、メッシュ図を示します。なおケース5と6(流速100 km/h、200 km/h)の計算では、境界層メッシュの第1層を厚くしています。
なお、本解析は乱流構造の調査が目的ではないため、抗力に影響が出ない範囲で、メッシュサイズを粗くして計算しています。メッシュサイズの影響については、補足②を参照してください。
-
図3．メッシュ図（メッシュ数：約52万）
4．計算結果
表2に、各ケースの抗力係数を示します。
表2．各計算ケースの抗力係数
補足：抗力係数は、抗力が安定した期間を平均して算出
【ケース1：流速 0.1 km/h】
【ケース2：流速 1 km/h】
*）第二不変量Q：回転と変形状況を示す指標で、値が大きいほど回転成分が卓越していることを示します
【ケース3：流速 10km/h】
【ケース4：流速 50 km/h】
【ケース5：流速 100 km/h】
【ケース6：流速 200 km/h】
図5に、球の抗力係数の実験データに、OpenFOAMの計算結果(赤丸)をプロットしたグラフを示します。
-
図5．実験データとの比較
参考資料：Prepared by Professor J. M. Cimbala, Penn State University Latest revision: 11 January 2012
5．まとめ
本解析事例では、オープンソースCFDソルバーの代表格であるOpenFOAM(Foundation版)を活用し、ベンチマークとして球体の空力特性および抗力係数の計算を実施しました。解析の結果、臨界レイノルズ数近傍での抗力係数の急激な低下（抗力危機：Drag Crisis）の現象を、実験データ（Experimental Data）の傾向を良好に再現することができました。
また、様々な商用のCFDソフトウェアと比較した場合、OpenFOAMの最大の強みはオープンソースソフトウェアである点が挙げられます。HPCやクラウドを用いた大規模な並列計算（Parallel Computing）であっても、コア数に応じたライセンス費用を一切気にすることなく、計算を実行することができます。
計算コストと解析精度を両立するOpenFOAMは、自動車、建築、ドローンなど、複雑な形状の空力設計や形状最適化（Optimization）において、商用ソフトと並ぶ有力な選択肢として、十分な実用性を持つ強力なソリューションとしてご活用いただけます。
補足②：メッシュサイズについて
球から剥離した乱れが遷移して渦構造を形成していく様子は、メッシュサイズに大きく依存します。本解析の目的は抗力係数の比較であるため、抗力に影響が出ない範囲でメッシュサイズを粗くして計算しましたが（総セル数：52万）、渦構造に着目した解析では、メッシュサイズを細かくする必要があります。参考として、総セル数を320万、1150万とした場合のアニメーションを示します。
【総セル数：52万】
【総セル数：320万】
【総セル数：1150万】
図6．メッシュサイズによる速度勾配テンソルの第二不変量の変化
補足③：速度勾配テンソルと第二不変量 Q
2次元の円柱の抗力計算で渦を可視化する際は、例えば、XY平面上で計算した場合、Z軸周りの渦度で描画すると、円柱から剥離する渦の様子を確認できます。
一方、3次元である球から剥離した乱れが遷移して形成される渦構造は、回転軸が様々な方向を向いており、時々刻々変化します。このため、渦度や流線で可視化してしまうと、渦構造が変化していく様子が捉えにくくなります。3次元の渦を描画する方法の一つとして、速度勾配テンソルの利用が挙げられます。
速度勾配テンソルは、流速ベクトルU の位置ベクトルx の勾配です。X,Y,Z軸方向の流速成分u、v、wに対してx、y、z方向の勾配を取るので、速度勾配は、3×3のテンソルとして表現されます。
速度勾配テンソル
この速度勾配テンソルには、固有値の特性によって、座標系を変更しても値が変化しない3つの不変量が存在します。それぞれ以下に示す流体の局所的な特性を示します。[3]
第一不変量I：流体の圧縮状態を示します(流体の収縮・発散)。非圧縮性流体の場合、0になります
第二不変量Q：回転と変形状況を示します
第三不変量R：湧き出しと吸い込み状況を示します
3次元の渦の描画では、第二不変量Qを利用します。速度勾配テンソルを、対称成分である速度テンソルSと非対称成分である回転テンソルWに分けて整理すると以下になります。
$$ Q=\frac{1}{2}\left(W^2-S^2\right) $$
第二不変量Qは、その値によって、WとSのどちらが支配的かを示していると見なせます(⇒Qが正の領域は回転が卓越)。第二不変量Qを、見やすくなる特定の正の値の等値面(Iso-Surface)で描画すると、球体から乱れが剥離し、渦構造を形成する様子を確認することができます。
なお、OpenFOAMでは、systemフォルダ下のcontrolDictファイル内に定義するfunctionオブジェクトで、Qを結果データとして追加出力することができます。
また、Paraviewでは、以下の手順でQの等値面を描画することができます。
1. Filtersメニュー ＞ Alphabetical > Contourをクリック
→Pipeline BrowserにContour1が追加されます
2. Contour ByをQに指定
3. Value Rangeに正の値を入力
4. Coloringで、等値面に着色したい結果項目を指定(補足：単色で着色したい場合はSolid Color)
5. Applyボタンをクリック
