「テーマ」
・Openfoamのlaminar(層流モデル)及びRAS(レイノルズ平均乱流モデル)計算における境界層の妥当性の確認
「方針」
OpenFoamを用いて境界層解析を実施し、粘性流体力学の理論式との比較を行ってみます。
尚Openfoamを用いた2D翼の抵抗係数確認にてpisoFoam+RAS(k-epsilon)よりもsimpleFoam+RAS(k-OmegaSST)の方が風洞試験値と一致することを確認しています。それを受けて本稿でもsimpleFoam+laminar及びsimpleFoam+RAS(k-OmegaSST)にて計算しています。
１．層流境界層
１．１ 理論式
層流境界層の挙動を把握するため、ここではBlasius方程式(下式)を選択します(参考文献[1]7.4項)。
・∂ux/∂x+∂uy/∂y=0 (質量保存の式)
・ux∂ux/∂x+uy∂uy/∂y=ν∂2ux/∂y2 (ナビエストークスの式）
ここで以下の変数変換を行います。
・η=y/√(νx/U)：境界層y方向(高さ)の無次元量
・f(η)=ψ/√(νUx)
但しψは流れ関数で以下を満足します。
・∂ψ/∂y=ux
・∂ψ/∂x=-uy
結果、ナビエストークスの式が以下の式に変形されます。
2f”‘+ff”=0：Blasius方程式
初期条件：f(0)=f'(0)=0、f'(∞)=1
本微分方程式は様様な解放があると考えますが、ここではScilabを使って解いてみました。使用したモデルとf'(η)の計算結果を下図に示します。理論式の境界層厚さδ99=4.91√(νx/U)(参考文献[1]7.5項)との比較よりScilab計算結果は妥当と考えます。
尚初期値f”(0)が未定のため、f'(∞)=1となるような初期値を探す必要があるのですが、今回は参考文献から借りてf”(0)=0.332057(参考文献[1]7.4項)としました。 f'(η)=ux/Uという関係式よりuxを求め、次項openfoamの計算結果と比較しています。
１．２ 解析モデル＆OpenFoam設定
層流境界層の解析を行うためのモデルを設定します。
・空間メッシュ
4.5(長さ)×0.6(高さ)×0.1(幅)ｍの直方体を270×40×1の合計10800マスに分割。
高さ方向は上下で20倍の長さのグラジュエーション。
・動粘度：空気 1.5e-5[m2/s]
・速度分布：一定流2[m/s]及び10[m/s](代表長さ２mでレイノルズ数2.7e+5及び1.3e+6に対応)
・境界条件：底面0～0.5[m]=slip、底面0.5～4.5[m]=noSlipで厚み0の薄板を表現
幅方向はtype emptyにより2D計算
・解析ソルバー：SimpleFoam(非圧縮・定常乱流解析ソルバー)
・乱流モデル：laminar(層流)。
・fvSchemes及びfvSolutionは以下を使用
\tutorials\incompressible\simpleFoam\motorBike
１．３ OpenFoam計算(層流境界層)
U=2[m/s]における計算結果を以下に示します。
境界層は非常に薄いので全体が見えるように横軸スケールを1/50にしました(スケールと速度表示も1/50になっています)。更に薄板端から0.1、2、3.5[m]の3点での縦方向速度分布及びBlausis方程式の数値解をエクセル図で示します。 結果、oenfoamの計算結果とBlasiusの数値解が非常によく一致しています。
また薄板に働く粘性摩擦力も確認してみます。
OpenFoamの計算結果より薄板の各要素面に働くせん断応力(wallShearStress)を面積積分することにより摩擦抗力係数CDfricを求めました。
CDfric=(ΣAi×τi)/(0.5U2*A)=0.00192
一方Blausis方程式の数値解より(参考文献[1]7.6項)
CDfric=1.328/√Re=1.328/(2[m/s]*4[m]/1.5e-5)=0.00182
結果、双方ほぼ同じ値を示しています。
尚本計算のy+は[Min3.40 Max12.1 Ave4.53]となっています。壁関数を使わないため、この程度としました。
上記計算は平板中央で局所的レイノルズ数Re＝2.7e+5の計算としています。このレイノルズ数は円柱、球の幾つかの実験値で層流から乱流に切り替わるRe=3e+5に合わせています。実空間ならば薄板前半で層流境界層、後半から乱流境界層に切り替わることが期待できる設定にしました。
そこで更に高いレイノルズ数で乱流を確実に発生させることを目的に速度10[m/s](Re＝1.3e+6)を計算してみました。結果を以下に示します(1/50スケールに圧縮しています)。
こちらもoenfoamの計算結果とBlasiusの数値解が非常によく一致しています(速度が増した分、境界層は薄くなっています)。また摩擦抗力係数CDfricもOpenfoam=0.000835、数値解=0.000813とほぼ一致しました。しかし期待したような乱流への遷移は起きていません(層流設定なので当然ですが)。
以上よりOpenfoamのlaminar設定では実験で起きる層流から乱流への遷移は再現出来ないということが分かりました。
尚本計算のy+は[Min2.86 Max13.1 Ave3.79]となっています。
また速度が5倍になった分、境界層が薄く、それに対応するため高さメッシュを2倍の80マス、グラジュエーションを50倍にしています。
付け足し１：参考に垂直方向の速度分布と流線を下図に示しておきます。境界層開始エリアで垂直方向に速度成分が発生しますが渦になるわけではなく、層流として後方に流れていく様子が見て取れます。
付け足し２：本論とはずれますが、風速2[m/s]と10[m/s]で抵抗係数が0.00182から0.000813と約1/2に減っています。これは感覚とずれますので要注意と感じます(表面上の速度勾配は大きくなり、抵抗そのものは増えますが、比較対象となる動圧が速度の2乗で増えています)。
２．乱流境界層
２．１ 理論式
海では上空から見ると大きな波が観察されますが、少し高度を下げるとそれよりも小さな波が観測されるそうで、これが何回も繰り返されるそうです。乱流境界層内部の渦も同じ構造のようで、大小の渦を全て把握することは大規模計算機でも容易ではないそうです。
その対応の一つとして速度uを平均速度Uとそこから偏差u’に分けて、平均速度の挙動を計算する方程式をRANS(Reynolds-averaged Navie-Stokesequations)と呼ぶそうです(参考文献[2]4.7.2項)。
しかしこの方程式は層流境界層の時のように机上計算で解くことは難しいそうです。
ここでは代わりに以下の経験式をOpenfoamの検証用に用います(参考文献[2]4.8.1項及び[3]13.6項]。
・境界層厚さδ
δ=0.37x*(Rex)-1/5=0.37x*(Ux/ν)-1/5
・境界層内速度分布
u/U=(y/δ)1/7
・摩擦抗力係数(c:全長4[m])
CDfric=0.1184*(Rec)-1/5=0.1184*(Uc/ν)-1/5
２．２ 解析モデル＆OpenFoam設定
乱流境界層の解析を行うためのモデルを設定します。
層流境界層のモデルとの違いは乱流モデルのみです。
・空間メッシュ
4.5(長さ)×0.6(高さ)×0.1(幅)ｍの直方体を270×40×1の合計10800マスに分割。
高さ方向は上下で20倍の長さのグラジュエーション。
・動粘度：空気 1.5e-5[m2/s]
・速度分布：一定流2[m/s]及び10[m/s](代表長さ2mでレイノルズ数2.7e+5及び1.3e+6に対応)
・境界条件：底面0～0.5[m]=slip、底面0.5～4.5[m]=noSlipで厚み0の薄板を表現
幅方向はtype emptyにより2D計算
・解析ソルバー：SimpleFoam(非圧縮・定常乱流解析ソルバー)
・乱流モデル：RAS k-OmegaSST。
・k、nut及びomegaファイルの設定は適当(速度2[m/s]及び10[m/s]&I＝1%で参考文献[4]4.5.5項に従う。)
・fvSchemes及びfvSolutionは以下を使用
\tutorials\incompressible\simpleFoam\motorBike
２．３ OpenFoam計算(乱流境界層)
U=2[m/s]における計算結果を以下に示します。こちらも全体が見えるように横軸スケールを1/50にしました(スケールと速度表示も1/50になっています)。
更に薄板端から0.1、2、3.5[m]の3点での縦方向速度分布及び経験式の数をエクセル図で示します。
結果、経験式との比較は傾向は概ね一致していますが、境界層厚さはOpenfoamの方が薄いようです。
一方摩擦抗力係数CDfricはOpenfoam=0.00494、経験式=0.00847と2倍弱の違いが出ています。これはOpenfoamの速度分布が経験式よりもy=0付近の勾配が緩いため摩擦せん断力が減ったためと推定されます。
尚本計算のy+は[Min6.86 Max12.2 Ave7.65]となっています。ここでは壁関数に頼っていないのでこの位でも問題ないと考えます。
念のため速度10[m/s]についても計算しておきます。以下に計算結果を示します。
結果、境界層厚さは2[m/s]の時よりも薄くなっています。経験式との比較は傾向は概ね一致していますが、境界層厚さはこちらもOpenfoamの方が薄いようです。
また摩擦抗力係数CDfricはOpenfoam=0.0036、経験式=0.00614と２倍弱の違い(2[m/s]の時とほぼ同じ)が出ています。
尚本計算のy+は[Min7.63 Max13.2 Ave8.33]となっています。
またこちらも速度が5倍になった分、境界層が薄く、それに対応するため高さメッシュを2倍の80マス、グラジュエーションを50倍にしました。
現状、Openfoamのk-OmegaSSTは経験式と十分な一致を示しているとは言い難いところです。そこで乱流モデルの基本であるk-epsilonを適用してみたのですが、結果は摩擦抗力係数CDfricはk-epsilon=0.0041とk-OmegaSST=0.0036とあまり変わりませんでした。
以上乱流境界層の計算結果をまとめてみます。
・概ね一致するが以下の違いが見られた。
・境界層厚さは経験式に対しOpenFoamが数割小さい
・摩擦抵抗も経験式に対しOpenFoamが約半分(危険側)層流境界層と違い、乱流境界層は簡単に試験値と一致するものではない様子が伺え、試験値とのコリレーションが重要になると考えます。
付け足し１：参考にｋ分布と流線を下図に示しておきます。境界層開始エリアで垂直方向に速度成分が発生します。一方境界層内では平均速度は渦になるわけではなく、層流のごとく後方に流れていく様子が見て取れます。但しｋ分布を見ると変動成分としての渦を認識出来ます(注意：k成分は渦度とせん断力項の二つが含まれるためｋ分布の全てが渦度を示している訳ではありません)。
付け足し2：本論と外れますが、壁関数はもともと粗いメッシュでも近似的に境界層を模擬する目的があります。しかし2.3項のメッシュはかなり細かく、実用上計算が重くなります。そこでここではy+が30～300位になるようにメッシュを敢えて粗くし（1/8倍）、その状態でどのような計算結果になるかを観察してみます。
速度10[m/s]における計算結果を示します。本計算のy+は[Min13.3Max68.9Ave56.6]となっており、壁関数としては適切な数値と考えます。
摩擦抗力係数CDfricは0.00329、細かいメッシュモデルが0.0036なので十分な一致と考えます。
結果、メッシュが粗いわりに目的である摩擦係数はメッシュが細かい時とそれ程変わらない結果が得られました。これは計算時間を考えると非常に便利と考えます(simpleFoamは早いのでそれ程のご利益はないのですがpisoFoamでは大きく効きます)。
３．層流境界層+乱流境界層
３．１ 解析モデル＆OpenFoam設定
Openfoamのk-OmegaSSTLMモデルは層流境界層と乱流境界層を同時に扱うそうなので試しに使用してみます。
乱流境界層のモデルとの違いは乱流モデルとそれに伴う初期ファイル設定のみです。
・空間メッシュ
4.5(長さ)×0.6(高さ)×0.1(幅)ｍの直方体を270×80×1の合計10800マスに分割。
高さ方向は上下で50倍の長さのグラジュエーション。
・動粘度：空気 1.5e-5[m2/s]
・速度分布：一定流10[m/s](代表長さ2mでレイノルズ数1.3e+6に対応)
・境界条件：底面0～0.5[m]=slip、底面0.5～4.5[m]=noSlipで厚み0の薄板を表現
幅方向はtype emptyにより2D計算
・解析ソルバー：SimpleFoam(非圧縮・定常乱流解析ソルバー)
・乱流モデル：RAS k-OmegaSSTLM(γ-Reθモデル)
γ(間欠係数)とReθ(運動量厚さレイノルズ数)よりモデル遷移を計算する
・k、nut及びomegaファイルの設定は適当(10[m/s]＆I=1%で参考文献[4]4.5.5項に従う。)
・gammaIntは全域0.1(≒層流)を設定、ReThetatはOpenFoamのHPヘルプに従う(下式)。
・fvSchemes及びfvSolutionは以下を使用
\tutorials\incompressible\simpleFoam\T3A
・1.3項及び2.3項と違い計算が収束しない場合がある。その場合は変数がおよそ一定になった段階で終了とした。
３．２ OpenFoam計算(層流+乱流境界層)
U=10[m/s]における計算結果を以下に示します(スケール1/50)。
U_yとｋ(乱流運動エネルギー)を観察することにより前半は層流で後半に乱流になる様子が見て取れます(注意：ｋ分布の全てが渦を示している訳でありません)。
また乱流開始地点はおよそx=2[m](Y=2.5[m]から端部距離0.5[m]を引く)となっており、その点でのレイノルズ数はRex=(Ux/ν)=(10*2/1.5e-5)=1.3e+6となり、一般的に乱流になるといわれているレイノルズ数1e+5～1e+6と一致しています。
また摩擦抗力係数CDfricはk-OmegaSSTLM=0.00231とk-OmegaSST=0.0036よりも小さくなっています。これは層流部分の抵抗が小さいからと考えられます。
尚層流から乱流に遷移するメカニズムは複数の要因が推定されます。したがって風洞試験における遷移点の位置も様様な因子によりばらつくと予想されます。対してk-OmegaSSTLMモデルでは初期設定ReThetatを適当に選ぶことで遷移点が前後することを確認しました。具体的には上図計算ではｋ=3/2*(U*1%)^2を用いていましたが、これを5%にした場合を以下に示します。
遷移点はおよそx=0.5[m]に前進しています。その点でのレイノルズ数はRex=(Ux/ν)=(10*0.5/1.5e-5)=3e+5となっています。初期条件を変えることにより乱流開始のきっかけの位置がかわってくるのだと思います。
４．まとめ
・simpleFoam+Laminarの組み合わせによりBlasiusの理論式と十分一致した。
・simpleFoam+Laminarでは乱流に遷移すべき速度でも層流を維持している。
・simpleFoam+RAS(k-ωSST)の組み合わせによる乱流境界層厚さは経験式よりも薄い。
このため摩擦抵抗係数も経験式よりも小さくなっている。
・simpleFoam+RAS(k-ωSST)では壁関数のおかげでメッシュを粗く(1/8倍)しても計算結果がそれほど変わらない。
・simpleFoam+RAS(k-ωSSTLM)の組み合わせにより層流⇒乱流の遷移が表現されている。
・simpleFoam+RAS(k-ωSSTLM)ではReθ(運動量厚さレイノルズ数)の初期値により遷移点が移動する。
Openfoamによる境界層解析はソルバーに強く依存することを理解しました。このため試験値とのコリレーションが重要となります。
参考文献
[1]流体力学 非圧縮性流体の流れ学 中山司著 森北出版
[2]航空宇宙工学テキストシリーズ 粘性流体力学 日本航空宇宙学会編 丸善出版
[3]航空宇宙工学テキストシリーズ 空気力学入門 日本航空宇宙学会編 丸善出版[4]OpenFoamによる熱移動と流れの数値解析 第2版 OpenFoamCAE学会編 森北出版
改訂記録
A改訂 2026.10.5 全面改訂(ソルバーをpisoFoamからsimpleFoamに変更)
