今回から円柱流れを考えてみる。
これまでのパターンではMapleで作成した解説を先に公表していたが、今回は他アプリケーションの結果をベンチマークとして、FlexPDEのFEM結果を比較してみたい。
1 ベンチマークとして利用したアプリケーション
FeatflowのHPを開いてBenchmarksをクリックすると、Single Phase Navier-Stokes Equationsのタイトルで円柱流れの結果が4件見つけられるはず。一番初歩的な「RE20, Laminar」のケースを今回は採用した。計算条件は以下の通り
（1）計算モデル形状
X-Y方向の2Dモデル。幅0.42m（Y）、長さL 2.2m （X）のダクトに流体を流す（深さ方向Zは1mとする）。ダクト流入口から0.2mの位置に、流れに抵抗するように半径r 50mmの円柱をZ方向に立てる。
（2）流体条件
密度ρ 1.0 kg/m3の気体。温度は規定なし。動粘度ν 0.001 m2/s（ρ=1.0なので、粘度μは0.001 Pa・sとなる）。ダクト入口から一様流速u 0.2 m/sで流し込む。
（3）境界条件
vx=0：ダクト壁、円柱表面
vy=0：ダクト壁、円柱表面
p=0：ダクト出口
2 FlexPDEの結果（FEM計算）
（1）モデルのグリッド（自動グリッドでノード105、セル165）...図1
流体は左側から流入する。
（2）修正レイノルズ数（MRe）...図2
カラーバーを見た通り全体的に約32.17なので、まだ層流状態である。ちなみにFeatflowベンチマークではRe20。
（3）流速コンター（vm、つまりsqrt(vx^2+vy^2)）...図3
円柱と壁の間が流速最大0.36m/sで、円柱後部は流速最低0m/sである事が判る。またダクト出口では0.28m/s最大で放物形に分布していることが図8で判明する。
先のベンチマークの流速コンター図で読み取れた値を下に表記する。概略近い値を示している。
FlexPDE Featflow
円柱ー壁間最大部（m/s） 0.36 0.4
円柱後部（m/s） 0 0
ダクト出口最大部（m/s） 0.28 0.28
（4）流速ベクトル（vx,vy）の拡大...図4
円柱後部で渦を巻いている、つまり双子渦の発生がはっきり判る。流体力学のテキストではレイノルズ数 6<Re<40の範囲で双子渦が見られると記載ある。
（5）圧力コンター（p）...図5
ダクト流入口の圧力は約0.07Pa、円柱先端部は0.13Paで最大圧力を示す。ダクト出口がゼロ（大気圧）になっているのは境界条件で指定しているからだ。
先のベンチマークの圧力コンター図で読み取れた値を下に表記する。概略近い値を示している。
FlexPDE Featflow
円柱先端部（Pa） 0.13 0.12
（6）圧力コンター（p）の拡大...図6
円柱左側に圧力の最大点0.13Paが、円柱最大径を過ぎた直後に圧力の最低点-0.01Paが見える。先のベンチマークの圧力コンター図も同様な絵である。
そのベンチマークには、円柱前後の差圧が0.1Pa実測されたことも記載されている。
ともかく圧力による抗力発生が容易に予想される絵だ。
（7）入口流速（vx）...図7
境界条件で一様流速0.2m/sを指示した通りになっている（ム？ギザギザじゃないの？と訝った方がいるかもしれないが、無料使用の範囲では精度がこの程度である）。
（8）出口流速（vx）...図8
0.28m/sを頂点に放物形の分布を示している事が判る。
（9）質量保存則（div_vのコンター）...図9
連続の式を計算した。カラーバーから、円柱付近を除き全体的にゼロと見て取れる。うまく発散して質量保存則が成立している証拠だ。
計算に使用したFlexPDEのコードとその説明は次回とする。
