「目的」
Opnefoamの回転チュートリアルを用いて遠心ポンプの基礎的流体解析を行う。
(インペラ回転がポンプ吐出流につながっていく様子を視覚で確認する。)
「方針」
一般的にはエネルギー効率論で済ましてしまうポンプ性能について、openfoamを用いた流体解析を行いて、視覚的に明確にしたいと考えます。
１．解析モデル
１．１ モデル方針
以下の形状を検討対象とします。
・半径125[ｍｍ]×厚さ10[ｍｍ]の容器
・半径25[ｍｍ]のハブ(中央部)
・75[ｍｍ]×10[ｍｍ]の放射状インペラ４枚
・流体：水。ハブから1[㍑/s]で吸入(予旋回なし)。
・インペラ回転1800[rpm]
１．２ OpenFoam設定
openfoamの設定を以下に示します。
・メッシュ
125[mm](半径方向)×360°(円周回り)x10[mm](高さ)を40x96x1の合計3840マスに分割。
インペラ周りは更に2倍に細分化。
またポンプ効率を計算するための検査面を半径75[mm]で設定(巻末補足3項)。
・積分幅：Δt=50μ[sec](初期値)
・クーラン数制限 0.5(計算開始時v*Δt/Δx=0.6x5e-5[m]/0.0017[ｍ]=0.02)
・回転速度1800[rpm](反時計回り)。
・ポンプ吸引量：1e-3[m3/s]。(type fixedValue; value nonuniform List<vector>)
・速度境界：インペラ表面はnonslip、外面はzeroGradient、上下面はempty条件。
・動粘度：水 8.57e-7[m2/s]
・解析ソルバーはSRFPimpleFoam(非圧縮・非定常乱流解析ソルバー)
・乱流モデルはRAS(kEpsilon)。
・k及びnutファイルの設定は適当(速度8[m/s]で参考文献[1]4.5.5項に従う。)
上記以外のパラメータは以下のOpenFoamチュートリアルの物をそのまま使用。 tutorials\incompressible\SRFPimpleFoam\rotor2D
本ソルバーはシングル・レファレンス・フレーム(SRF)機能を持ったpimpleFoamです。このソルバーは計算領域を回転させた状態を回転座標系から観測した状態で模擬します。具体的には計算領域全域にコリオリ力を印加します(巻末補足1項にてSRFの注意すべき特性を示します)。
尚本来インペラ回転により発生する負圧(サクション)により吸入量を決定するべきなのですが、このモデルでは初期条件ファイルにて固定しています(今回は1e-3[m3/s])。本来は吸入に関する流体回路をモデル化する必要があるのですが今回は大きく簡略化(吸入量固定)しています。
２．計算結果
２．１ ４枚インペラ遠心ポンプ
以下に回転座標系での速度Urelと圧力pの計算結果を示します。
インペラによって流体が撹拌されているようです。ただインペラを脱した後に流速が上がり、且つ圧力が上昇しています。この点はSRFPimpleFoamの問題と考えられ、巻末補足1項にて確認します。
もう少し流れの様子を確認するために流線を表示してみます。
なんとインペラによって加速された流れの一部が再びインペラ側に戻ってきてしまっています。具体的数値として吐出量1.6[㍑/s]のうち、0.6[㍑]が再度ポンプに戻ってきていました。吐出量の4割は利用出来ない流れとなっていました。なお総量としては吸引量1[㍑/s]と一致します。
検査面の流量の様子を切り出してみました。中央辺りから逆流している様子が分かります。
これらの図は以下のポンプの本質を示していると考えます。
・コリオリ力により吸入孔から入ってきた流体は後方インペラに積極的に向かっていく
・インペラによって流体に速度＆圧力を与える
・前方インペラ背面の負圧にコリオリ偏向が重なり、ポンプ外部からの逆流となる。
それではエネルギー効率を確認してみます。
インペラ直後の流量値から流体の得た仕事量Lout[W]を以下の式(参考文献[2]8.2項)から計算することが定番ですが、CFD計算結果は流れが均一でないため計算値の選定が困難となります。
L=Tω=ρQ(uvcosα)
そこで吸入孔とインペラの外縁に検査面を作成し、流体エネルギーの差を求めることにしました(検査面設定：巻末補足3項)。
結果は以下となります(それぞれ絶対座標系での面積平均を計算しています※１)。
・流速(＠流量)：
吸入孔 0.69[m/s]＠1e-3[m3/s]
検査面(吐出部) 12.5[m/s]@1.59e-3[m3/s]
検査面(逆流部) 12.2[m/s]@0.6e-3[m3/s]
・圧力
吸入孔 -165.6[Pa]：見易さからこの圧力を基準点に改めます。よって0[Pa]
検査面(吐出部) -79.1[Pa]：基準値変更により86.5[Pa]
検査面(逆流部) -79.3[Pa]：基準値変更により86.3[Pa]
流体の持つエネルギーはベルヌーイの式より
∴Lout=1/2*ρQ(Vo-Vi)+ρQ(Po-Pi)=165.2[W]
内訳 吸入孔 ：-0.5*1e3*1.0e-3*0.692-1e3*1.0e-3*0=-0.2[W]
検査面(吐出部) ：0.5*1e3*1.59e-3*12.52+1e3*1.59e-3*86.5=261.8[W]
検査面(逆流部) ：-0.5*1e3*0.6e-3*12.22-1e3*0.6e-3*86.3=-96.4[W]
（密度ρ：1e3[kg/m3]、流量Q：1e-3[m3/s])
続いてインペラが流体に与えた仕事量Lin[W]を計算します。
この計算にはopenfoamの汎用関数「forceCoeffs」のトルク係数CmPitchを用いました。CmPitchはインペラが流体に与えたトルクを示しています(詳しくは巻末4項に示します)。
openfoamにおけるforceCoeffsの計算結果が下図となります。このトルクT(グラフ収束値2.84e-3)に回転数ω(1800[rpm]=60π[rad/s])をかけることで仕事入力量Lin[W]を求めました。
Lin=1/2*ρ*CmPitch*ω＝0.5*1000*2.84e-3*60π=267.7[W]
上述効率値の妥当性を実験経験式(ウイスナーの式)により確認してみました。結果、61.3％となり上記値と十分な一致を見ました(巻末補足2項)。よってSRFPimpleFoamを用いた今回の計算は妥当と考えて良さそうです。
※１：回転座標系では流量保存が維持されていないようでした。一方絶対座標系では維持されていることを確認しました。そこで本計算では絶対座標系の速度を使用しました。
２．２ 流線形８枚インペラ遠心ポンプ
前項の検討より遠心ポンプの性能を向上させるにはインペラ部の逆流を防ぐことが重要と考えます。
そこでインペラ形状を流線形(後向き羽根)にし、かつ枚数を増やしてインペラ空間を小さくすることで逆流を減らしてみたいと考えます。
以下の形状を作成しました。これはインペラのない状態で600[rpm]における流体の流れを模擬したものです。但し補足１項に示す通り、本ソルバーではインペラのない状態の計算は問題を抱えています。しかしここでは厳密なインペラ形状の必要はないのでこのまま使用します。
前項と同様の計算を実施しました。結果が以下です。
・流速(＠流量)：
吸入孔 1.1[m/s]＠1e-3[m3/s]
検査面(吐出部) 12.2[m/s]@1.47e-3[m3/s]
検査面(逆流部) 12.5[m/s]@0.55e-3[m3/s]
・圧力
吸入孔 -149.6[Pa]：基準変更により0[Pa]
検査面(吐出部) -52.4[Pa]：基準値変更により97.2[Pa]
検査面(逆流部) -52.9[Pa]：基準値変更により96.7[Pa]
流体の持つエネルギーはベルヌーイの式より
∴Lout=1/2*ρQ(Vo-Vi)+ρQ(Po-Pi)=155.5[W]
内訳 吸入孔 ：-0.5*1e3*1.0e-3*1.12-1e3*1.0e-3*0=-0.6[W]
検査面(吐出部) ：0.5*1e3*1.47e-3*12.22+1e3*1.47e-3*97.2=252.3[W]
検査面(逆流部) ：-0.5*1e3*0.55e-3*12.52-1e3*0.55e-3*96.7=-96.2[W]
・仕事入力量Lin[W]：
Lin=1/2*ρ*CmPitch*ω＝0.5*1000*2.40e-3*60π=226.2[W]
結果ポンプ効率η=155.5/226.2=69%と放射状インペラの効率62％よりも上がりました。
しかし流線形8枚インペラの効率は思ったほどよくない結果でした。その原因として依然逆流の多さ(３割強)があると考えます。上図流線の様子を見ると依然前方インペラの背面に粗密領域が出来てしまっています。
そこで試しに回転数を1800[rpm]から半分の900[rpm]に変更してみます。
計算結果が以下です。
・流速(＠流量)：
吸入孔 1.1[m/s]＠1e-3[m3/s]
検査面(吐出部) 5.9[m/s]@1.2e-3[m3/s]
検査面(逆流部) 6.6[m/s]@0.2e-3[m3/s]
・圧力
吸入孔 -34.1[Pa]：基準変更により0[Pa]
検査面(吐出部) -9.6[Pa]：基準値変更により24.5[Pa]
検査面(逆流部) -10.1[Pa]：基準値変更により24.0[Pa]
流体の持つエネルギーはベルヌーイの式より
∴Lout=1/2*ρQ(Vo-Vi)+ρQ(Po-Pi)=40.5[W]
内訳 吸入孔 ：-0.5*1e3*1.0e-3*1.12-1e3*1.0e-3*0=-0.6[W]
検査面(吐出部) ：0.5*1e3*1.2e-3*5.92+1e3*1.2e-3*24.5=50.3[W]
検査面(逆流部) ：-0.5*1e3*0.2e-3*6.62-1e3*0.2e-3*24.0=-9.2[W]
・仕事入力量Lin[W]：
Lin=1/2*ρ*CmPitch*ω＝0.5*1000*9.70e-4*30π=45.7[W]
結果ポンプ効率η=40.5/45.7=89%と効率が十分上昇しました。
(但し吐出能力の絶対値(Lout)は1800[rpm]に対して約1/4になっているので嬉しくはないのですが) 。
3．まとめ
上述流体解析結果より遠心ポンプの働きは以下にまとめられると考えます。
・コリオリ力により吸入孔から入ってきた流体は後方インペラに積極的に向かっていく
・インペラによって流体に速度＆圧力を与える
・吸入孔から流入してきた流体は均一にインペラ空間に広がるのではなく、前方インペラ背面の負圧とコリオリ偏向の２つの要因により後方インペラ側に偏る。
・前方インペラ背面近傍は周囲に対し相対的に負圧になりポンプ外部からの逆流につながる。
上述が遠心ポンプの効率に大きく寄与していると考えます。
またこれらの動きから吸入孔にて、コリオリ力を弱める方向に流体に事前旋回を与えることでより効率が上がることが予想されます(予旋回と呼ばれるそうです)。
尚実際のポンプでは更に以下に示す損失を考えていく必要があると考えます。
・圧力漏れ(インペラ周辺での漏れ+インペラ前後のつり合い穴によるもの)
・流管回路損失(管路抵抗及び形状損失)
・渦損失(流体内渦の熱への変換分)
・機械的損失(ポンプベアリング周りの摩擦損失)
・原動機損失(電動モータ効率)
「補足１」SRFモデルの計算確認
まずはSRFモデルの問題点を示します。
中央ハブから絶対座標系で放射状に噴出した場合を絶対座標系で見た流れを下図に示します。
例え1800[rpm]で高速回転していても、インペラがないため絶対座標系では放射状に流れていかなければならないのに絶対座標系で見た流れが下図のように螺旋運動を示しています。
また本来コリオリ力は見かけの力のため物体の進路変更のみで増速には寄与しないはずです。しかし下図の計算結果は中央からの流れが螺旋運動しながら、増速しています。
このようにSRFPimpleFoamソルバーではインペラ等障害物のない空間では実世界では起きないコリオリ力による増速が起きてしまっているようです。
このSRFモデルの問題の原因を推定してみます(当事業所独自の推定という点をご留意下さい)。
推定原因
・コリオリ力を印加する際に有限要素における離散化により角度誤差が発生し、コリオリ力の一部が速度方向に分力として印加されてしまった。
この仮定のもと、コリオリ力の角度誤差で螺旋運動につながるかを大雑把な机上計算で確認してみます。
回転座標系での流体の流れを下図に示します。
・大雑把な速度変化の状況
中心部(ハブ) 円周方向 1800[rpm]×内側半径25[mm]=4.7[m/s]
周縁部(外周) 円周方向 5.6[m/s](図中より)
回転時間 4回転 周縁部で代表計算 4×.1(外側半径)[m]×2π/5.6[m/s]=0.45[sec]
∴加速度=(5.6-4.7)/0.45=2[m/s2]
・コリオリ力による加速度
周縁部で代表 2×ω×v=2111[m/s2]
ω：1800[rpm](60π[rad/s])、v：5.6[m/s]
結果、僅かtan-1(2/2111)=0.05°の角度誤差がopenfoam計算内で発生するとコリオリ力の分力により螺旋運動の説明が可能となります。これほどわずかな角度誤差は有限要素法では避けられそうにありません。
続いて圧力について確認します。
SRFモデルではモデル最周縁に構造物を設定していないのですが、周縁面より吐出可能な流量は吸入量(1e-3[m3/s])と同じとなり大幅に吐出量が制限されます。この吐出制限が仮想的な壁を演じます。このためSRFモデルの最周縁部では円運動を維持するための向心力(遠心力)が発生します。
この様子は宇宙ステーションにおける人工重力と同じで、重力ポテンシャルを形成します。
このポテンシャルを机上計算してみます。
p=1/(2πR*h)*∫mrω2=1/(2πR*h)*∫(ρ*2πrhdr)*rω2=ρω2/R*∫r2dr=58k[Pa]
ρ：水の密度1000[kg/m3]
ω：絶対座標系でのUy=13.3[m/s](図中)より1270[rpm]相当(133[rad/s])。
(回転数1800[rpm]で回転座標系は回っているけれど、流体はそこまで回転していない状態)
r：半径0.025～0.1[m]
R：外半径 0.1[m]
h：高さ0.01[m]
openfoam計算結果である中央部負圧-58k[Pa]とよく一致します。圧力に関しては速度と違い、実際の物理現象との相違はなさそうです。
尚モデル上は周縁面をp=0としたため、遠心力による圧力は吸入孔での負の圧力(-58k[Pa])として現れます(openfoamでの圧力はp/ρであり、密度分1000kg/m3をかける必要があります)。
以上からSRFモデルにおける注意点をまとめておきます。
〇速度に関し
・回転数が大きい場合、コリオリ力が相当大きくなる。
・FEM離散化に伴う角度誤差によりコリオリ力が流れの絶対速度を加速させてしまう。
・結果螺旋運動となる
よってSRFモデルを使用する際は流れが螺旋運動に入らないような構造物(インペラ等)を設定した方が良いと考えられます。
〇圧力に関し
・実際の物理現象との相違はなさそう。
・SRFモデルでは周縁の仮想的な壁により円運動に制限される。
・実際のポンプでもケーシングが存在し、そこで流れは円運動に制限される。
「補足２」経験式を用いた検証
本文2.1項にて計算されたSRFpimplefoamの効率62％の妥当性を下記経験式(参考文献[3]4.3項)により確認してみます。
Hth=1/g*(u2vu2)=1/g*u2{(1-k2)*u2-vm2cotβe}
Hth：理論揚程
k2：滑り係数
ウイスナー(Wiesner)の式より
k2=1-{1-1/z0.7*√(sinβe)}*{1-[((r1/r2)-ε)/(1-ε)]3}=1-.6211*.9872=0.387
ε=1/exp[8.16*(sinβe)/z] =0.13
z：インペラ枚数 4枚
r1：インペラ内側半径 0.025[m]
r2：インペラ外側半径 0.075[m]
βe：インペラ出口角 放射状のため90[deg] ⇒ cotβe=0
u2：円周速度 1800[rpm]×0.075[m]=14.14[m/s]
vm2：吐出速度(法線方向) 1e-3[m3/s]/(2π*0.075*0.01[m2])=0.212[m/s]
理論揚程HthにρQgを掛けると仕事量Lになります(ρ=1e3[kg/m3]、Q=1e-3[m3/s])。
またk2=0の時、効率100％(L0と表記します)となります。
以上より上記経験式から求められる効率は
L0=ρQ*u2^2＝200[W]
L=ρQ*u2^2*(1-k2)=122.6[W]
η=L/L0=61.3%
結果SRFPimpleFoam計算結果62％と十分な一致をみることが出来ました(ここまでの一致はたまたまだと思います)。
「補足３」検査面の設定
本文中ではポンプの出力効率を計算するため、インペラ径に合わせた円柱面を設定しました。この検査面の設定がチュートリアルには見当たらないため、ここにまとめておきます。
(但し当事業所はこのソフトウェアの専門家ではないので、その点をご留意下さい。)
尚検査面を設定すると計算結果に若干の変化が見られますが、無視出来る程度です。
Step1：形状作成
何らかのCAD(FreeCAD等)を用いて円柱面のSTLファイルを用意します。ここでは「cylinder.stl」という名前にしておきます。そのファイルをopenfoamの/constant/triSurfaceフォルダに格納します。
Step2：snappyHexMeshDict&createPatchDict＆controlDictの追加文
snappyHexMeshDict内に以下のようにcylinder.stlを定義し、コマンド「snappyHexMesh -overwrite」を実行します。これによりmeshに検査面が出てきます(paraview「xx.foam」propertiesの「Read zone」欄をチェックのこと)。
geometry{
cylinder.stl{
type triSurfaceMesh;
name insPect;} } <=検査面をinsPectと命名
refinementSurfaces{
insPect{
level (0 0);
faceType baffle; ＜＝表裏のある境界面を設定。
faceZone insPectZone;
patchInfo{type patch;}}} ＜＝流体の行き来を許可
次に後処理で検査面の性質をコマンド「createPatch -overwrite」で定義します。
尚本文モデルは2Dモデルにするため、このコマンドの前に「extrudeMesh」を実施しています。
pointSync true;
patches({
{ name insPect; ＜＝検査面表側の名前
patchInfo{
type cyclic; ＜＝master-slave面の宣言
neighbourPatch insPect_slave;}
constructFrom patches;
patches (insPect);}
{ name insPect_slave; ＜＝検査面裏側の名前
patchInfo{
type cyclic;
neighbourPatch insPect;}
constructFrom patches;
patches (insPect_slave); } );
更に初期条件ファイル0/Urel、p、ｋ、epsilon、nutに以下の検査面初期条件(全共通)を追加します。
insPect {
type cyclic;
neighbourPatch insPect_slave; }
inspect_slave{
type cyclic;
neighbourPatch inspect; }
また計算結果を記録するため「controlDict」に以下のfunctionを追加して準備が終了します。
functions {
insPectFaces{
type surfaces;
libs (“libsampling.so”);
executeControl runTime;
executeInterval 0.1;
writeControl runTime;
writeInterval 0.1;
surfaceFormat vtk;
fields(
p
U
phi
Sf );
surfaces (
insPectSurface {
type faceZone;
zones ( insPectZone );
interpolate false; } ); } }
以上の事前準備後にSRFPimpleFoamを実施すると以下のフォルダに検査面での計算値がvtpファイルとして保存されます。
/postProcessing/insPectFaces/1.0/insPectSurface.vtp
Step3：データ読込
paraviewにて上述vtpファイルを読み込みます。
これで目的のデータは入手出来たのですが、面平均等の後処理を実施するためにCSVデータを作成してみます。
以下の順にparaviewのコマンドを実施していきます。
・「Filters」>「Alphabetical」＞「Cell size」＞「Apply」：各面情報を追加
・「Filters」>「Alphabetical」＞「Surface Normals」＞「Apply」(Compute Cell Normals欄をチェック)：各面の面ベクトルの追加
・「Filters」>「Alphabetical」＞「Integrate Variables」＞「Apply」:各面データ欄がポップアップしてきます。
・データ欄のShowingを「Surface Normals1」に、Attribute「Cell Data」に選択することで各面の諸データ(p,phi,U)が下図のごとく表示されます。
・最後に右端のアイコン「Export spreadsheet」をクリックしてCSVファイルとして保存出来ます。
「追記１」
任意の面からの流量phiは以下の式で求めることが出来ます。
phi = A*(nx*Ux + ny*Uy + nz*Uz)
A：任意面の面積
nx、ny、nz：面の法線ベクトル
Ux、Uy、Uz：面における流速ベクトル
openfoamにおいては以下の変数を用いることで任意の面の流量を計算出来ます。
nx、ny、nz：Normals_0、Normals_1、Normals_2
Ux、Uy、Uz：U_0、U_1、U_2
「追記２」
回転座標系と固定座標系の変換では接線方向の速度のみ反転する必要があります。
この場合、任意の面で行う変換は流速ベクトル(Ux、Uy)を時計回りに面の傾きθ度回転し、そこでy成分を反転、再度反時計回りにθ度回転することで得られます。ベクトル行列で表現すると以下になります。
但し2次元のみ適用の式となります(下式は行列計算のつもりです)。
|Ux’| = |cosθ -sinθ| |1 0| |cosθ sinθ|| Ux |
|Uy’| = |sinθ cosθ| |0 -1||-sinθ cosθ|| Uy |
(Ux‘、Uy‘)：反転後流速ベクトル
openfoamにおいては以下の変数を用いることで任意の面での接線方向流速のみ反転することが出来ます。
(本文ではこの式を用いて回転座標系での初期速度を計算しています。)
cosθ、sinθ：Normals_0、Normals_1
Ux、Uy：U_0、U_1
「追記３」
step3のデータ読込はZoneでなくてもモデル中の任意の面の計算結果を抽出出来ます。よって次項補足3で示します汎用関数「forceCoeffs」を用いなくても対象物に発生する荷重の抽出が可能となるのですが、扱うデータは大量となります。
補足４：汎用関数「forceCoeffs」の利用
インペラが行った仕事量Lin[W]の計算にはopenfoamの汎用関数「forceCoeffs」を用いました。
この関数を実行すると以下の変数の時間変化が得られます。
・揚力係数 Cl
・抵抗係数 Cd
・横力係数 Cs
・ピッチ、ロール、ヨーモーメント係数CmPitch、CmRoll、CmYaw
この関数は使用前に以下のパラメータを設定します。
rhoInf 1000; 密度
liftDir (0 1 0); 揚力方向
dragDir (1 0 0); 抵抗方向
CofR (0 0 0); 軸原点
pitchAxis (0 0 1); ピッチ軸
magUInf 1.0; 速度絶対値
lRef 1.0; 代表長さ
Aref 1.0; 代表面積
そこでこの関数のパラメータを以下とすることでCmPitchからインペラにかかるトルクTを計算しました(これは当事業所のオリジナル処理という点をご留意下さい)。
・magUInf=1
・IRef=1
・Aref=1
・rhoInf：適当(計算中で使われていない模様)
T=1/2*ρmagUInf2*CmPitch*Aref*IRef
∴T=1/2*CmPitch
参考文献
[1]OpenFOAMによる熱移動と流れの数値解析第２版 一般社団法人オープンCAE学会編 森北出版
[2]改訂 流体の基礎と応用 盛田泰司著 東京電機大学出版局
[3]SI版 渦巻ポンプの設計 設計製図の基礎 柏原 俊規 他共著 パワー社
HA-technodefy 