月刊下水道 2015 11 月号 VOL.38 No.13

有限要素法(FreeFem++)による三次元流体解析－ 手軽に流れを観察する その２ －中日本建設コンサルタント（株）

中根進

1. まえがき
筆者は、本誌 Vol.36.No.10 (2013 年増刊号)で「格子ボルツマン法による下水流れの可視化－手軽に流れを観察する－」と題して、フリーソフト(Blender：Fluid)を使って自由水面の大きな変動を扱える三次元格子ボルツマン法とその解析事例を紹介した。
このソフトでは、自由水面の変動を可視化、数値出力はできたものの、筆者の技量の範囲で流体内部の流れを可視化できなかった。
筆者が、文献１）を参考に三次元に拡張したエクセルベースの差分法による三次元流体解
析は、あまりにも計算時間を要することから、本稿では、下水道施設の流体内部の流れを流速ベクトルや流線で可視化することを目的に有限要素法のフリーソフト「Fleefem++」を使い三次元流体を解析することとした。
本稿では、このフリーソフト「FreeFem++ Ver.3.34」および３次元解析領域の生成と有限要素メッシュを切るために使用したフリーソフト「Gmsh Ver.2.9.1」を紹介し、計算結果を FreeFem++に実装されている「medit とフリーソフト「Paraview Ver.4.3.1」で可視化した下水道施設における解析事例を示す。
2. 流体基礎式（１）ナビエ・ストークス方程式
FreeFem++を使って有限要素法による３次元流体解析を始めるにあたり、流体の流れを表すナビエ･ストークスの方程式を示す。
ナビエ・ストークスの方程式（Navier-Stokes equations）１）は、非圧縮性粘性流体の運動を記述する２階非線形偏微分方程式である。
連続の式(密度一定)

ｘ方向の運動方程式

非定常項

粘性(拡散)項

移流項

圧力項

非線形の効果慣性力：流体が慣性によってこれまでの運動をそのまま続けようとする力

外力項

粘性力：
液体のもつ粘性が運動を均一化する力
圧力勾配による力：
圧力変化があることによって生ずる力1ｙ方向の運動方程式、ｚ方向の運動方程式省略(文献１）参照)

ここに、ｕ：ｘ方向の流速Ｆｘ：単位質量あたりに働く外力(重力、表面張力など)
ρ：密度ｐ：圧力

μ 1
ν＝ ＝
ρ Ｒe

ν：動粘性係数μ：粘性係数
Ｒe：レイノルズ数（２）ナビエ・ストークス方程式の弱形式

有限要素法では、ナビエ・ストークス方程式の偏微分方程式を弱形式といわれる方程式にする。FreeFem++には、その弱形式の方程式に対応する関数をもっている。FreeFem++上でその関数と境界条件を加えてプログラムし、方程式を解き、各有限要素の流速、圧力を得る。
本 稿 で は 、 FreeFem++ の ホ ー ム ペ ー ジ に 例 題※

解析領域(三次元形状)の作成※
計算

面番号の設定時間発展

グラフィック
画面確認誤差(examples++-3d)としてアップされ、ナビエ・ストークス方程式を弱形式化してある NSI3d.edp を利用して解析した。
解析領域のメッシュ分割

結果確認

グラフィック画面確認

ファイル出力3. FreeFem++による流体解析
3.1 流体解析の手順
FreeFem++による解析の手順は、図-１の通りである。
Gmsh を使い三次元解析領域の定義とメッシュの作成まで行う。(図-１の一点鎖線の範囲)
非線形偏微分方程式に対する数値計算をFleeFem++により行い、代数方程式の解（数値解）を求め、 時間発展に対しては差分法を使用する。

メッシュファイルの読込

ParaViewによる可視化

有限要素空間、空間内変数の定義
偏微分方程式の定義

計算結果は、FreeFem++に実装されている medit や別のフリーソフトである Paraview によって可視化する。

境界条件の設定※

この有限要素メッシュに対して、FreeFem++の中で有限要素空間や空間内の流速や圧力など変数を定義する。
2

図-１ FreeFem++による解析手順

方程式を解くために、流入面、流出面、水面、壁面などの境界条件や時間ｔ＝0 の初期条件を設定する。
3.2 Gmsh による三次元メッシュの作成
FreeFem++における三次元領域の作成には、文献２）の３つの方法と直感的で操作しやすいGmsh を使う方法がある。
ⅰ) TetGen を使う方法ⅱ) コマンド Buildlayers を使う方法ⅲ) 三次元の有限要素データのファイルをコマンド readmesh3 で読み込む
Gmsh は、入力用のダイアログから点の座標値を入力し、点と点の間を結んで(マウスでクリックして)線とする。線を左廻りで結んで面を作り、複数の面を選択して三次元領域を作成する。
3.3 Paraview による解析結果の可視化
FreeFem++には実装された medit を使い解析結果を３次元で可視化することができる。
しかし、FreeFem++で流線を描くためには、流れ関数を導入して、流速ベクトルとは別に計算しなければならない。
Paraview には Stream Tracer というフィルタが用意されていて、流速ベクトルさえあれば、流線を描くことができる。上述したように三次元ではメッシュ数が多いと計算できないので、筆者は今のところ、粗いメッシュで計算して、Paraview で補間計算して流線を描画するようにしている。
FreeFem++には、計算した流速ベクトル値を savevtk()というコマンドを使って拡張子.vtk のファイルに出力できる。Paraview でこのファイルを読み込み、Paraview の機能を使って流線などを描く。
4．数値計算条件
4.1 境界条件の設定
偏微分方程式で表された流体の方程式は、初期条件、境界条件を与えて解を求める。不適切な境界条件では計算できない場合や、現実離れした結果になる場合もある。
流体の方程式を解く場合には適切な境界条件を設定する必要がある。
方程式は時間微分を含むため、ある指定された時間(たとえばｔ＝0)における境界条件も必要となる。この境界条件を特に初期条件１）という。
１)水面、壁面の境界条件
壁面や水面の境界では流体がその境界から離れることはできないので、その表面上で、表面の法線方向の相対的な速度が 0 にする。
壁面境界条件：x 方向の流速ｕ＝0，ｙ方向のｖ＝0，ｚ方向のｗ＝0２)流入面、流出面の初期条件
3

方程式を解くために、流入面、流出面のどちらかに流速(ｕ,ｖ,ｗ)などの条件を与える。
流入面、流出面両方に流速などの条件を与える方程式を解くことができなくなる。
4.2 クーラン条件
時間発展でナビエ・ストークスの方程式を解くので、時間刻みΔｔとメッシュの大きさが適切でないとエラーがでて計算が止まる。メッシュの大きさと時間刻みΔｔの比をクーラン数 C といい、１以下にする必要がある。
C ＝ｕ × Δt/Δx, Δｙ, Δｚ ≦ １
ここで、
C：クーラン数
ｕ,ｖ,ｗ：流速Δt：時間刻み
Δx,Δｙ,Δｚ：格子の大きさ5．数値計算事例
5.1 ポンプ井解析例
ポンプ井とポンプ吸込部をモデル化し、ポンプ井の壁やポンプ周りの流れを可視化し、ポンプに空気や真空域の吸込がないよう、渦や負圧になりやすい位置を推定する。
Gmsh で作成した流体領域を図-２に示し、メッシュを切った状態を図-３に示す。

図-２ ポンプ井解析領域の作成(Gmshによる表示)

図-３ 三次元メッシュ分割例(Gmsh による表示)

Gmsh で 作 成 し た 図 - ３ の メ ッ シ ュ をFreeFem++に読み込んで、medit で表示すると図-４となる。
Gmsh では図-４のメッシュよりも、細かくすることができるが、FreeFem++で実行すると、方程式を解く際に使用するソルバーがメモリーオバーとなり、計算ができなくなる。
4

図-４ ポンプ井メッシュ(Freefem++の描画ツール medit による表示)

壁面・水面・ポンプ吸込口の境界条件、初期条件を以下のように設定する。
壁面 ：ｕ＝0、ｖ＝0、ｗ＝0
水面 ：ｗ＝0
流入面： 初期条件としてポンプ吸込側で流速を与えるので、流入面は、Neumann(ノイマン)境界条件３）として、流速を与えない。
ポンプ吸込口：
ポンプや機械撹拌では、羽根の回転により流れには傾きがあるため、ディリクレ(Dirichlet)境界条件３）として以下の計算式で与える。
ポンプ吸込部周りの流速、流向の設定を図-５に示す。

図-５ ポンプ吸込部の流速、流向の設定(Graph-R Ver.2.32)による表示

計算モデルの各部寸法を文献４）に示す値で設定した。口径 800mmφの各部寸法のポンプ井に対して、1000mmφのポンプを想定し、ポンプ井の壁やポンプ周りに有害な真空域(流速が著しく小さい領域)や渦などが表現できることを期待した。
FreeFem++の描画ツール medit を使い、流体表面の流速分布を図-６に示す。

図-６ポンプ井の流速分布結果(Freefem++の描画ツール medit による表示)

Paraview を使って、流体内部の流線を表示すると図-７、８となる。
5

図-７ ポンプ井の流線(Paraview による表示)

図-８ ポンプ吸込部の流 線 (Paraviewによる表示)

5.2 機械撹拌解析例
ポンプ吸込部と同じプログラムで、吸込部を回転撹拌部に変え、ｚ方向の流速ｗをポン
プ吸込と反対側に与えた。水面を時計回りに機械撹拌することを想定した。入力した撹拌
機周りの流速、流向の設定を図-９に示す。

図-９撹拌機の流速、流向の設定(Graph-R)

撹拌槽モデルの大きさは、幅２ m×長４ m×水深３ m、表面撹拌径 1.0 m であり、その概要を図-10、メッシュを図-11 に示す。
撹拌

図-11 機械(表面)撹拌メッシュ

図-10 機械(表面)撹拌モデル

Paraview を使って、流線を表示すると図-12、13 となる。図-13 は、流線方向を矢印で表示した。
6

図-12 機械(表面)撹拌による流線

図-13 機械(表面)撹拌による縦断面の流線とその方向各要素の流速を抽出し、流速の頻度分布を作成する。活性汚泥法の反応槽では、活性汚泥フロックの沈殿防止のため、底部流速を 0.1 m/sec 以上を確保するとされている。汚水、雨水の沈砂池では、砂分を沈降さ5せる目的で池内流速を 0.3 m/sec

4.5以下としている。
(図-14)から判断すると、0.1 m/sec以下のゾーンが 39.4％ある。槽の底部に分布しているかどうかは、底部の流速を出力して調べる必要があるが、0.2m/sec 以下も 72.8％

39.4％

3.53頻度

上記撹拌槽は、流速の頻度分布

0.1 m/sec

40.2 m/sec

2.5

72.8％

21.5
10.5000.20.40.6 0.8
11.2 1.4
各部の撹拌速度 m/sec

1.61.82

図-14 槽内の流速の頻度分布あり、撹拌強度がこの例では不足しているように思われる。
5.3 曝気槽解析例

汚水の流入、流出を持ち、空気撹拌による槽内の流速分布を解析する。空気撹拌は、空気の上昇にともなう気泡と液体の摩擦によって流れが発生するもので、底で発生した気泡は、水面近くで水圧が小さくなると気泡径が大きくなり摩擦力も変化すると考えられる。数値解析モデルとしては、気液混相流として解析する必要があるが、ここでは、散気面に圧力差Δｐを与えて流れが発生する散気面ようにした。
図-15 全面曝気方式解析モデル

7(１)全面曝気方式
解析モデルを図-15 に示す。
槽寸法：幅 6.0 m×長 10.0 m×水深 5.0 m
散気面：幅 3.0 m×長 9.0 m
槽内の流線は、図-16 となり、撹拌されている様子がうかがえる。図-16 左の流線は、槽の中のある線上(ｙ＝0(槽中央)、ｚ＝4.0(水深))を通過する流線を表している。槽内流速の頻度分布を図-17 に示すが、0.1 m/sec 以下となる領域が 26％程度あった。流速が小さい領域は、鉛直方向、水平方向の断面の流速を出力して探す必要がある。

図-16 全面曝気による槽内の流線
14旋回流式,開口中央
12

確率密度

1086420.340.26000.10.2
槽内流速 m/sec

0.30.4

図-18 旋回流の槽と散気面のメッシュ

図-17 槽内の流速の頻度分布(２)旋回流方式
槽寸法：幅 6.0 m×長 10.0 m×水深 5.0 m
散気面：幅 1.2 m×長 9.0 m
槽と散気面のメッシュを図-18 に示す。槽内の流線を図-19、20 に示す。図-20 は流線に流速の大きさを着色したものである。
槽内流速の頻度分布は、図-24 であり、全面曝気と同程度の流速が確保されている。
図-19 旋回流の流線(paraView による表示)8

14全面曝気,開口中央
12

確率密度

1086420.27000.360.10.2
槽内流速 m/sec

0.30.4

図-21 槽内の流速の頻度分布

図-20 旋回流の槽横断面の流線

5.4 分水槽解析例
汚水を均等に分配させることを目的とした図-22、23 の分水槽を解析する。既報の格子ボルツマン法では、堰からの流出量を把握するには、図-24 に示す数値モデルの中に計量槽を作り、この中の水面座標を出力する必要があった。そこで、プログラムの容易な有限要素法で解析し、分水量を算定したものである。

1.0

図-27 に示す。
各ステップ水路の可動堰からの流出量(分配)割合を表-1 示す。その結果、均等な割合で分配され0.61.20.90.25
1.20.25

初沈流出睡0.6し、分配割合を算定した。越流部の流速ベクトルを1.62.75

0.25 3.60

図-22 分水槽平面図ていることが明らかになった。

格子ボルツマン法

計量槽図-27 格子ボルツマン法による分水状況
図-26 分水槽解析モデル

91.2

可動堰越流部の流速をテキストファイルで出力0.43.8線を図-29 に示す。

導流壁

共通水路

場合の流体表面の流速分布を図-28 に示す。また流1.5から、共通水路から v＝0.108 m/sec で流入させた

初沈流出水0.6

図-23 に示す初沈流出水を水面上(v＝0.4 m/sec)

2.6

図-25 分水槽の流速分布(medit による表示)
図-26 分水槽内の流線(paraView による表示)

表-１ 分水槽流出口の分配割合
ステップ水路名 流量 －分配割合

ステップ１11.20.33

ステップ２11.50.34

ステップ３11.30.33

計

34.01.00

図-27 分水槽流出口の流速分布(Graph-R による流速)

6．まとめ
筆者が FreeFem++で３次元流体解析ができるようになった手順を示しているもので、解析結果から具体的に問題点を洗い出し、検討したものではない。既報「－手軽に流れを観察する－」の第２弾として、三次元の有限要素法を紹介したものである。また、Paraview による流線の矢印による表現、流線を速度により着色した表現など多彩な表現の一部を紹介した。
下水処理施設での流体解析例として、ポンプ井、反応槽、撹拌槽、分水槽などの流況状況を示すことができた。ポンプ井、撹拌槽では、水の流れとして流線を数値解析で示し、流速分布から止水域と思われる割合や場所を示すことができた。分水槽の解析では、分配割合を示すことができた。これら数例でも皆さん自身が持つ下水処理施設における水理的な問題を解決できるのではないかと思われる。
なお、下水処理施設での流体は、非定常な流れであり、水処理施設にあっては乱流、密度流や気液混相流として扱う必要がある。また汚泥処理施設の汚泥では、非ニュートン流体として扱うなど、本稿で紹介した数値解析は利用しても意味のないものかもしれない。
市販ソフトであれば、これらの下水の特性にあった流れを解析できるので、購入するか、フリーのソフトで開発して、問題の解決にあたっていただきたい。また、紹介した計算のパラメータはレイノルズ数(動粘性係数)だけなので、実施設での調査や水理模型実験などで検証することが必要である。
10

＜参考文献＞
１）河村哲也著：エクセル・シミュレーション入門 山海堂 2004 年 4 月 22 日 pp.51～56
２）日本応用数理学会監修、大塚厚二･高石武史著：有限要素法で学ぶ現象と数理－FreeFem++数理思考プログラム－ p.83３）樫山和男（中央大学）：有限要素法による流れ解析の基礎（２）-Navier-Stokes 流れhttp://www30u-toyama.ac.jp/okumura/fem/pdftext/first-day/Kashiyama.pdf
４）一般社団法人 河川ポンプ施設技術協会：揚排水ポンプ設備技術基準・同解説
年２月 ｐ.3-2311平成 27