第 11 章 非粘性流れの FEM 解析
非粘性流れの数値解析は，コンピュータ以前に手計算で始められている．その代表的手法は緩和法，特性の方法，渦法である．緩和法 (relaxation methods) は，Laplace 方程式，Poisson 方程式などの境界値問題を差分法で解く際に得られる大型疎行列の連立 1 次方程式を，効率よく解く方法である．Southwell の緩和法の教科書には，1 冊めにこれらの方程式の解法，2 巻めに重調和方程式やその他の方程式の解法が述べられている．特性の方法 (methods of characteristics) は，非定常１次元流れや２次元超音速流れなどの双曲型
微分方程式の初期値問題を図式解法で解くもので，衝撃波や滑り面を正確に計算することができる．特性曲線法 (characteristic curve methods) とホドグラフ特性法 (hodograph characteristic methods) がよく知られている．また渦法はポテンシャル流れ問題の解法で，一様流中にある物体を渦分布で置き換えその流れを解析するものである．特にパネル法 (panel methods) はコンピュータの時代になって飛行機の設計に用いられ，機体表面と後流中に渦を分布させ機体まわりの流れを解析するものである．
コンピュータの時代に入って非粘性流れ問題は，差分法，有限要素法などの各種の数値解法で解かれるようになった．非粘性流れの基礎方程式は，基本的には連続方程式と Euler の運動方程式で，圧縮性流れの場合にはこれにエネルギー方程式が加わる．ただし２次元流れは，定常また非定常でも非圧縮性流れの場合には流れ関数を定義することができ，その 2 階楕円型微分方程式として定式化することができる．また渦なしのポテンシャル流れは，速度ポテンシャルの Laplace 方程式 (非圧縮性流れ) または Poisson 方程式 (圧縮性流れ) として定式化することができる．Laplace 方程式で支配される流れは境界積分法により境界のみを対象に計算することもできる．
この章では２つの非粘性流れ問題を異なるアプローチにより今日広く使われている有限要素法 (FEM) で解くことにする．すなわち始めに定常２次元翼列流れが，流れ関数を導入し，その楕円型微分方程式の境界値問題を解くことによって求められる．次に３次元亜音速翼列流れが，速度ポテンシャルの Poisson 方程式の境界値問題を解くことによって求められる．

1

2

非粘性流れの FEM 解析

11.1

楕円型微分方程式

FEM の解法に入る前に本節では，Laplace 方程式や Poisson 方程式などの楕円型偏微分方程式 (elliptic
partial differential equations) の解の性質について簡単に述べる．
微分方程式は，独立変数の数，従属変数の数 n，階数 k ，線形性，型により解の性質が異なり，解法も違ってくる．n 個の k 階方程式は一般に kn 個の 1 階の方程式の系に置換えることができる．ここでは詳細は省略するが1 ，この方程式系の特性方程式の根がすべて共役複素根の場合に，その方程式ないし方程式系
は楕円型 (elliptic) といわれる．楕円型微分方程式にはポテンシャル問題の Laplace 方程式，Poisson 方程式，弾性問題に登場する重調和方程式などがある．線形非線形の定義はやや曖昧である．linear, semi-linear,quasi-linear までを線形とし，Monge-Ampère 方程式のような genuine non-linear のみを非線形とするものもある．線形方程式は重ね合わせの原理を適用して解くことができる．また半線形と準線形方程式は，相当の線形方程式を被摂動方程式に摂動法で解くことができるので，重ね合わせの原理は広義の線形方程式に適用できることになる．数値解法では，重ね合わせの原理や摂動法という言葉はほとんど出てこないが，このような考えは反復計算の過程で多用されている．
k 階の線形楕円型微分方程式 L(u) = 0 を考える．楕円型方程式は通常 境界値問題として解かれ，閉じた境界上に，境界値として関数自身とその k−1 階までの法線微分 u, ∂u/∂n, . . . , ∂ k−1 u/∂nk−1 の中から k/2
個のものが与えられる．ただし ∂/∂n = n ·∇ は法線 (方向への偏) 微分，n は境界上の外向き単位法線ベク
トルである．例えば 2 階の Laplace 方程式の境界条件は，すべての境界上の 関数 u またはその法線微分

un = ∂u/∂n の値である．しかしこれはあくまでも基本を述べたもので k/2 ずつ与えない場合もある．更には無限遠方や周期条件を含む多連結領域への対処法も必要である．
本節では主にポテンシャル問題について述べる．Laplace 方程式と Poisson 方程式は次のように表わされる．
4u = 0(11.1)

4u = −ωn µ(x)(11.2)

ただし 4 ≡ ∇2 は Laplace 演算子で，デカルト座標系では２次元で ∂ 2/∂x2 + ∂ 2/∂y 2 ，３次元で ∂ 2/∂x2 +∂ 2/∂y 2 +∂ 2/∂z 2 となる．また ωn は立体角で２次元で ω2 = 2π ，３次元で ω3 = 4π となる．微分方程式の解
は，幾何学的に xu 空間内の積分曲面で表わされる．この曲面の高さは関数値 u，その傾きは関数の微分値
ux , . . . ，また曲率は 2 階微分を表わしている．Laplace 方程式の解すなわち調和関数 (harmonic function)
は至るところ平均曲率 0 の曲面で，また Poisson 方程式 4u = f の解は，f > 0 のところで平均曲率正，f < 0 のところで平均曲率負の曲面になる．
Laplace 方程式 (11.1) は，ある点 ξ から x までの距離 r = x−ξ のみによる解

1
1


log
(n = 2)
2πru(x) = γ(r) =(11.3)
1
1


(n
≥
3)

(n−2) ωn rn−2を持つ．この解は r = 0 にある特異点 (characteristic singularity) 周りのポテンシャル，物理的には質量 ωnの質点周りの重力ポテンシャル場を表わしている．Poisson 方程式 (11.2) に関しては，右辺の µ(x) が領域

1 16.1 節または 17.1 節参照

楕円型微分方程式

G で有界で積分可能ならば，
Z1log µ(ξ) dξru(x) = ZG1µ(ξ) dξ

G r

3

(２次元)

(11.4)(３次元)

で表わされる関数 u(x) とその 1 階微分は一様連続になる．更に µ(x) が G で区分的連続微分可能ならば

u(x) の 2 階微分も G で連続になり，この u(x) は Poisson 方程式 (11.2) の解になる (証略)．この解は 無限に広がる空間内の 領域 G における密度 µ(x) の質量分布によるポテンシャル場を表わしている．
グリーンの公式 (Green’s formulas)：

次図に示すように x 空間内の開領域を G，その境界を Γ ，その外向き単位法線ベクトルを n とする．Γ は区分的に滑らか，なお G + Γ は閉領域とする．Gauss の定理
RR∇X dx = Γ n X ds に X(x) = · v(x)∇u(x) を代入すれば次式が得られる．
GZ
Z
Z∇u·∇v dx + v 4u dx = vun ds(11.5)GGΓ

ただし un = ∂u/∂n = n ·∇u，u(x) は G で 2 回連続微分可能，G+Γ で連続微分可能，v(x) は G で連続微分可能，G+Γ で連続な関数とする．また同様の条件のもと次式が成立する．
ZZ
Z∇u·∇v dx + u 4v dx = uvn dsGGΓ

これらの式の差を取れば次のグリーンの公式が導かれる．
ZZ(u4v−v4u)dx = (uvn − vun )ds
G(11.6)

Γ

ただし u(x) と v(x) は領域 G で 2 回連続微分可能で領域 G+Γ で連続微分可能な関数とする．
式 (11.5) に 4u = 0, v = 1 を代入すれば Gauss の積分定理 (Gauss’ integral theorem) が得られる．
ZZ
un ds =
∇u·dss = 0(11.7)
ΓΓ

ただし dss = n ds は Γ 上の微小面ベクトルである．この式は，関数 u(x) が G で正則調和で G+Γ で連続微分可能ならば，上記の法線微分 un の積分が 0 になるという連続の条件を表わしている．
次に式 (11.5) に u = v, 4u = 0 を代入すれば Dirichlet 積分 (Dirichlet’s integral) が得られる．
ZZ(∇u)2 dx = uun dsGΓ

ここでも関数 u(x) は G で正則調和で，G+Γ で連続微分可能とする．このとき• Γ で u = 0 ならば u は G で恒等的に 0 になる．

(11.8)

4

非粘性流れの FEM 解析• Γ で un = 0 ならば u は G で一定になる．
Laplace 方程式の解曲面の性質を考えればがってんできよう．
ここで一息入れる．関数 u(x) が一様連続の ‘一様’ は x の定義領域の至るところという意味である．‘連続’の基本は Lipschitz 連続で，Lipschitz 連続は不等式 |u(x1 )−u(x0 )| < κ|x1 −x0 | によって判定される．x0 は¡¢任意の点，x1 はその近傍点で，x1 −x0 の大きさを十分小さくしたときに勾配 u(x1 )−u(x0 ) /(x1 −x0 ) の大√きさが十分大きい数 κ よりも小さければ関数 u は点 x0 で連続と判定される．この定義では関数 u(x) = x
は区間 x ≥ 0 で繋がっているが x = 0 の近傍で不連続と判定される．この不都合を忌避するため Hölder 連
続では不等式|u(x1 )−u(x0 )| < κ|x1 −x0 |α ,0<α≤1(11.9)

によって判定が行われる．関数 u(x) が式 (11.9) を満足するときに，u(x) は α 次の Hölder 連続であるとい√う．区間 x ≥ 0 で x は，0 < α ≤ 1/2 のときに式 (11.9) を満足するので 1/2 次の Hölder 連続である．区
分的連続 (piecewise continuous) 関数は，通常いくつかの不連続を持つ有界な連続関数である．

n 回連続微分可能 (n-time continuously differentiable) な関数は，n 階までの微分が存在しそれらがすべて連続になる関数である．領域 G で n 回連続微分可能な関数 u(x) は，u ∈ C n in G のように表わされる．
例えば領域 G で求めた解 u(x) がその境界 Γ 上の境界値 u = φ に接続する条件は u ∈ C 0 in G+Γ のように表わされる．無限回連続微分可能な関数を解析関数 (analytic function) という．これに対し少なくも当面必要な回数までは連続微分可能な関数を正則関数 (regular function) という．解析的であれば正則である
が，正則であることは一般に解析的ではない．またここでは滑らか (smooth) は解析的と同義とする．十分
滑らか (sufficiently smooth) は一般には更に滑らかなことをいうが，ここでは当面必要なだけは滑らかという意味とする．数学者は連続性をやかましくいうが，読者は上記の例からどの程度連続性をうたえばよいのか容易に推測できよう．他方 数値解法では連続性にはほとんど言及しないが，例えば解析関数すなわち

Taylor 級数の収束を前提とする差分式をところ構わず使うのは問題であろう．

平均値の定理 (mean value theorem)： 上図に示すような点 P を中心に持つ同心円または同心球の間の領
域 G を考え，グリーンの公式 (11.6) に 4u = 0, v = log(1/r) (２次元), = 1/r (３次元) を代入し，ガウスの積分定理 (11.7) を用いれば，左辺は 0 になり，次式が得られる．
ZZ
1
1
u ds = lim
u ds ≡ uP(２次元)
R0 →0 2πR0 Γ
2πR Γ
0
Z
Z
1
1
u
ds
=
lim
u ds ≡ uP(３次元)
2
R0 →0 4πR0
4πR2 Γ
Γ0
ただし R, R0 は円または球の半径である．これよりある点 P の調和関数の値 uP は，P を中心に持つ円または球の表面における u の値の算術平均に等しい．ただし関数 u はこの円または球内で正則でその閉包楕円型微分方程式

5(closure) G+Γ で連続とする．また点 P の調和関数の値 uP はこれらの円または球の内部の u の値の算術平均に等しいということもできる．
最大値の原理 (maximum principle)： 平均値の定理より，関数 u(x) が連続領域 G で正則調和で，その閉包で連続ならば，u の最大値と最小値は常に境界 Γ 上にある．u の最大値または最小値が u の内部にあるのは u が定数値のときのみである2 ．速度ポテンシャル φ に対しては，

u=0∇4φ = 4(∇φ) = − 4u
ただし u は流速である．これより ポテンシャル流れの流速の最大値と最小値は必ず境界上にあることがわ
かる．Poisson 方程式 4u = f の場合には，f > 0 ならば u(x) の最大値は常に境界 Γ 上にあるが最小値は必ずしも Γ 上にあるとはいえないことになり，また f < 0 ならばその最小値は常に境界 Γ 上にあるが最大値は必ずしも Γ 上にあるとはいえないことになる．

ある領域の内部，境界上，または外部の点 x を考える．今上図に示すようにその領域から点 x まわりの微小領域を除いた領域を G とし，この領域 G に対するグリーンの公式 (11.6) に v = log(1/r) (２次元),= 1/r (３次元) を代入すれば，次式が得られる3 ．
Z
I
I
1
1
1
1
1
∂ ³ 1´u(x) = −4u(ξ) log dξ +un (σ) log dσ −u(σ)
log dσ
(２次元)
p G
r
p Γ
r
p Γ
∂n
r
(11.10)
Z
Z
Z
1
1
1
1
1
∂ ³1´
u(x) = −4u(ξ) dξ+un (σ) dσ
−
u(σ)
dσ
(３次元)
p G
r
p Γ
r
p Γ
∂n r


(x inside G) ωnp = ωn /2(x on Γ )0
(x outside G)R
ただし r = ξ − x， Γ f (σ)dσ は G の外側境界上の積分，ωn は立体角で２次元で 2π ，３次元で 4π である4 ．
この式は，u ∈ C 2 in G，u ∈ C 1 in G+Γ ならば• 領域 G 内の密度 − 4 u/ωn = µ の空間分布• 境界 Γ 上の密度 un /ωn = ρ0 の single-layer
2 G の内部の点に最大値があるとする．その点を中心に G 内に円または球を取れば，u が定数値でなければ最大値はこの円または球面上にあることになる．この新たな最大値に対しても同様のことが言え，結局 最大値は境界上にあることになる．
3 ３次元の x
in の場合には，内側境界 Γ0 上の積分は
Z
Z
1
∂ “1”
lim
un (σ) dσ = 0, limu(σ)
dσ = 4πu(x)
r→0 Γ
r→0 Γ
r
∂n r
0
0
4 これらの式は本来，点 ξ に分布する質量等による点 x のポテンシャルを求めるもので式 (11.4) のように書く方が素直である．しかしこのように定義した r を用いる式も多く見られる．これらの式は，ξ にある質点による x のポテンシャルの値が x にある同じ質
量をもつ質点による ξ のポテンシャルの値に等しいこと，これらのポテンシャルの微分も大きさ等しいが負号が異なることを考えて作られたものと解釈される．

6

非粘性流れの FEM 解析• 境界 Γ 上の密度 −u/ωn = ρ の double-layer
によるポテンシャルと考えられる．µ は有界で積分可能，ρ, ρ0 も同様．single-layer はその線 (２次元) または面 (３次元) を横断して関数自身が連続でその法線微分値が不連続，また double-layer は関数が不連続でその法線微分値が同じになる線または面である．境界 Γ 上の点 0 に着目し，この点に接する領域 G 内の点を上添字 −，この点に接する領域外の点を上添字 + を付けて表わせば
0−u+
n − un = −ωn ρ0 ,• single-layer：

u+ = u− = u0 ,• double-layer：u+ − u0 = u0 − u− =

ωn
ρ0 ,
2

−
u+
t = ut = ut0−u+n = un ,−u+
t − ut0 = ut0 − ut =

ωn
ρt0
2

である．ただし下添字 t は線または面の接線方向の微分を示す．一般に Poisson 方程式の境界値問題では，
領域 G 内の方程式 4u = f と 境界 Γ 上の u または un の値が与えられる．したがって境界 Γ 上の u の値の与えられる部分では un の値，また un の値の与えられる部分では u の値が未知になり，式 (11.10) はそのままこの境界値問題の解を与えるものではない．

u が速度ポテンシャル φ の場合には，領域 G で 4φ = 0，また境界 Γ の速度が与えられれば境界における流れの法線速度 φn が既知になる．特に境界が静止している場合には式 (11.10) は次のようになる．
I1∂ ³ 1´φ(x) = −φ(σ)
log dσ
(２次元)
p Γ
∂n
r
Z
1
∂ ³1´
φ(x) = −
φ(σ)
dσ
(３次元)
p Γ
∂n r
この問題では Γ 上に double layer すなわち渦層を置くことによって流れ場を構成することができる．渦層の渦の強さは Γ 上の法線速度 φn がすべて 0 になるように決定される．また境界が動いているところでは，

φn が既知の法線速度になるように決定される．
ここで具体例として一様流中に置かれた翼型を過る２次元流れを考える．この翼型周りには通常循環流が存在し次式で定義される循環は一般に 0 にはならない．
II
C = u ·dss = − ∇φ·dss = φs − φtただしHは翼型を左回りに一巡する閉曲線に沿う積分で，添字 s は積分の起点，t は終点を意味しこれらは同一点である．この流れ場は２連結で，φ の境界値問題の解を一価関数で表わすためには，通常 翼面から無限遠まで切断 (cut) を入れ，ここでポテンシャルの値を跳躍 (jump) させることが必要である．翼面上の流れでは φn = 0 したがって翼の内側でも φn = 0，翼の内部では φ 一定になる．例えば翼の内部で φ = 0，
翼後縁から切断を入れ，後縁の切断の上側で φ = 0 と置く．その下側の φ の値すなわち切断を横切っての跳躍量は Kutta の条件，翼後縁から流れが滑らかに流れ去るという条件を満足するように調整される．な
お当然のことながら翼面上の流速 −φt は，境界上ではなく境界に接する流れの内部で求められる．
次に一様流中に置かれた翼周りの３次元流れを考える．この流れでは通常 翼端渦と後流 (随伴) 渦面が考慮される．これらの渦は翼を down wash するもので無視できない．後流渦面は翼端渦間に翼負荷がスパン方向に一様でないときに生じるもので，翼後縁から脱落する渦糸によって構成されるひとつの流面である．
この面の上下の流速は大きさ等しく方向の異なるものになる．切断は後流渦面に沿って入れられ，Kutta の条件もかされる．この流れは通常パネル法などの渦法によって解かれている．次に翼列を通る流れを考える．この流れにおいても翼 (羽根) の周りに循環が生じるので切断を入れる必要があるが，通常 陽にではな
く，翼列下流領域の周期条件にその役割を担わせている．また Kutta の条件は流出角の調整によって満足

楕円型微分方程式7

させることができる．３次元翼列流れにおいても後流渦面は考慮すべきである．渦法で２次元翼列流れが解かれている．

Laplace 方程式または Poisson 方程式の境界値問題は，構造格子を用い差分法や有限体積法で，また有限要素法では非構造格子を用いても解かれている．ディリクレ問題 (Dirichlet’s problem) は，境界条件 u = ϕ(s)on Γ を満足する微分方程式 4u = f (x) の G で正則な解 u(x) を求める問題である．ただし f (x) は区分的連続関数で，u ∈ C 0 in G+Γ とする．またノイマン問題 (Neumann’s problem) は境界条件 un = ϕ0 (s) on Γ ，u = ϕP at x = xP ∈ Γ を満足する微分方程式 4u = f (x) の G で正則な解 u を求める問題である．ただし
f (x) は区分的連続関数で，u ∈ C 1 in G+Γ とする．なお混合境界値問題 (mixed boundary value problem)
は，境界 Γ のある部分に関数値 u = ϕ(s)，残りの部分に法線微分値 un = ϕ0 (s) の与えられる上記と同様の問題である．
非圧縮性流体の定常２次元平面流れは，流れ関数 (stream function)ψ をZx;ψ = ψ0 + i z · u ×dx
ux = ψy , uy = −ψxx0のように定義し，流れ関数の Laplace 方程式または Poisson 方程式の境界値問題を解くことによって求めることができる．この流れ関数は渦あり流れや粘性流れにも適用できる．また非定常の非圧縮性流れ，定常の圧縮性流れ，更には流面が平面でなく例えば回転面 (revolutional surface) の場合，それに付加される流路高さの変化する場合にも定義できる．なお流れ関数を３次元に拡張したベクトルポテンシャルを用いる
アプローチも提案されているが，式の数多くソレノイダル条件の問題もあり，ほとんど用いられていない．
流れ関数の値は流線に沿って一定で，したがって静止境界上でも一定になる．また運動境界上の流れ関数の値はその定義にしたがって変化する．
一様流中に置かれた翼型を過る流れは，翼型から十分離れた遠方境界に一様流の流れ関数の値，翼型表面にある流れ関数の値 ψB を与え，流れ関数の方程式のディリクレ問題を解くことによって解析することができる．ψB の値は Kutta の条件を満足するよう予測子修正子法によって決定される．また翼列を通る流れ
は，一つの翼間流路とその上流と下流への延長部分からなる計算領域に対し，翼型表面にこの部分流路の流量に合わせて流れ関数の値を与え，上流と下流領域の上下の境界には流れ関数とその法線微分の周期条
件をかし，上流境界には流れの流入角，下流境界には流出角に合わせ法線微分値を与え，流れ関数方程式の混合境界値問題を解くことによって解析することができる．なお流出角はあらかじめ与えられるものではなく Kutta の条件を満足するように予測子修正子法によって調整される．詳しくは 11.3 節のプログラムとその説明参照．
次にターボ機械の回転している羽根車を通る３次元流れを考える．羽根車の相対流れは渦あり流れにな
るが，絶対流れはその上流にある静止流体が直接的に流れ込む場合を想定すれば渦なし流れになり，速度ポ
テンシャル φ を定義することができる．この場合の流れは φ の Laplace 方程式の周期条件を含む混合境界値問題を解くことによって解析することができる．羽根の下流側には単独翼の場合と同様に後流渦面を考慮すべきである．この境界値問題の計算領域は，一つの翼間流路とその上流と下流領域への延長部分からなり，特に下流領域の周期境界は後流渦面に合わせて取られる．また境界条件は，ハブ，ケーシングまたは羽u，ただし u は羽根車の周速，上流境界には φ
根車のシュラウド面上では φn = 0，羽根面上では φn = ± n ·uの値，下流境界には φn = −wm ，ただし wm は子午線流速，が与えられ，また周期境界には φ と φn の周期条件がかされる．なお後流渦面，下流境界の循環と子午線速度の分布は反復計算の過程で決定されるものである．その一つの計算例を 11.4.4 項に示す．
グリーン関数 (Green’s function)： グリーン関数を用いる解法は古典的数学の教科書には必ず登場するが普通 数値解法に用いられることはない．ある領域 G に対する 4u = 0 のグリーン関数 K(x, ξ) は，G 内8非粘性流れの FEM 解析の点 ξ に特異点を持つ固有基本解で次のように表わされる．

K(x, ξ) = γ(r) + w(11.11)

ただし γ(r) は式 (11.4) で与えられる x = ξ に特異点を持ちこの点を除いて 4γ = 0 になる関数，また w はG+Γ で連続，G で 4w = 0 で，Γ 上でグリーン関数を K(x, ξ) = 0 または Kn (x, ξ) = 0 にするための関数である．グリーン関数は対称性があり K(x, ξ) = K(ξ, x) である．なおこのラプラス方程式に対するグリーン関数の定義は一般の楕円型方程式に，また非同次の境界条件にも拡張されている．
グリーン関数の Laplace 方程式または Poisson 方程式の境界値問題への応用： 初めに Dirichlet 問題，境界条件 u = φ on Γ を満足する微分方程式 4u = f の G + Γ で連続で，G で正則な解 u を求める問題を考える．式 (11.6) に v = K を代入し，式 (11.7) を用いれば，
ZZ
Zu(x) = − K 4u dξ − Kn u dσ + Kun dσGΓΓ

更に 4u = f in G，u = φ on Γ ，K = 0 on Γ を代入すれば，この問題の解がグリーン関数を用い次のように表わされる．
ZZu(x) = − K(x, ξ)f (ξ) dξ − Kn (x, σ)φ(σ) dσG(11.12)

Γ

次に Neumann 問題を考える．境界条件を un = φ0 on Γ のように与え，Kn = 0 on Γ の方のグリーン関数
を用いれば，この問題の解は次のようになる．
ZZu(x) = − K(x, ξ)f (ξ) dξ + K(x, σ)φ0 (σ) dσ + const.
G(11.13)Γ

ただし const. は与えられた境界条件 uP = φ at P ∈ Γ から決定される．この解法の問題点は，一般的境界形状に対しグリーン関数を求めることが境界値問題の解を求めるに等しい労力を要することである．
ポアソンの積分公式 (Poisson’s integral formula)： ２次元の円または３次元の球の内部領域 G で 4u = 0，
その表面 Γ で u = φ なる Dirichlet 問題の解は，φ が連続ならば，次の Poisson の積分公式で与えられる．
Z
Rn−2 (R2 −ρ2 ) φu(x) =
dωn(11.14)ωnrn

r2 = ρ2 +R2 −2ρR cos θ
ただし R は下図に示すように円または球の半径，ρ はこれらの中心から点 x までの距離，r は点 x から Γ
上の点 ξ までの距離である．この公式は次のように導出される．円または球形の領域 G に対するグリーン関数は，点 P の特異点とその相反点 (reflection point)P1 の特異点の和で次のように与えられる．
³ ρr ´1
K = γ(r) + w = γ(r) − γR9

楕円型微分方程式

原点 O から相反点 P1 までの距離は R2/ρ，前図の 4OP Q と 4OP1 Q は相似で r : ρ = r1 : R になる．この関数 K は境界 Γ 上で 0 になりグリーン関数である．次にこの関数の Γ 上での法線微分を求めればKn = cos(∠OQP )∂∂ ³ρr1´ R2 −ρ2 0γ(r) − cos(∠OQP1 ) γ
=
γ (r)
∂r
∂r1 R
Rr

ここでは次の関係が用いられた．

R2 +r1 2 −(R2/ρ)2R2 +r2 −ρ2, cos(∠OQP1 ) =
,
2Rr
2Rr1
1
ρ
∂ ³ρr1´ρ ∂∂0γ(r) = −γγ(r) = γ 0 (r)
≡
γ
(r),
=
∂r
ωn rn−1
∂r1 R
R ∂r
R

cos(∠OQP ) =この Kn を式 (11.12) に用いればポアソンの積分公式が得られる．

Poisson’s integral formula•?
analyticity ofharmonic function•?

mean value theorem•?

maximum principle
•

?

??

uniqueness theorem??

?

Weierstrass’ convergence
theoremcompactness theorem??

reflection principle

existence theoremof solution

ポアソンの積分公式に由来する重要なもを以下に述べる (上図参照)．ポアソンの積分公式 (11.14) で点 Pを円または球の中心 O に取れば，すなわち ρ = 0, r = R と置けば調和関数の平均値の定理が得られる．また Laplace 方程式や Poisson 方程式のディリクレ問題が異なる解を持つとする．これらの解の差に最大値の原理を適用すれば異なる解は持ち得ないことになり，この問題の一意性の定理 (uniqueness theorem) が導かれる．
調和関数の解析性 (analyticity of harmonic function)：ある領域 G で調和で正則なすべての関数は Gの内点の近傍でべき級数に展開できるので，解析的である (証略)．この定理により，ポテンシャル流れの理論では，正則 (regular) と解析性 (analytical) が同義に扱われることがある．
反射の原理 (reflection priciple)： ある領域 G で調和で，その閉包 G+Γ で連続な関数 u は境界の球面または平面の部分 Γ0 で u = 0 ならば，u はこの部分の境界を越えて反射によって解析接続することができる (証略)．
以上述べた Laplace 方程式や Poisson 方程式で記述されるポテンシャル問題に関する定理は次のように一般化される．一般の 2 階の楕円型微分方程式は次のように表わされる．

L[u] ≡nXi,k=1

aik uxi xk +

n
X
i=1

bi uxi + cu = f(11.15)10

非粘性流れの FEM 解析

ただし領域 G 内で f, aik , bi , c ∈ Cm+α とする．このとき 2 階連続微分可能な解 u は G 内で u ∈ Cm+2+αである．更にその滑らかな境界 Γ 上の境界値 φ ∈ Cm+2+α を満足する解 u は G+Γ で u ∈ Cm+2+α である．２次元の Laplace 方程式の場合には，a11 = a22 = 1, a12 = b1 = b2 = c = f = 0 であるから，すべての係数 ∈ C∞ になり，解 u は解析的ということになる．また例えば２次元 Poisson 方程式の右辺 f が領域

G を 2 分する境界線の右側で 0，左側で 1 の場合には，この境界線上で f ∈ C−1 になり，解 u ∈ C1 ，すな
わち u とその 1 階微分は連続で，2 階微分は不連続ということになる．この場合には境界線の右側で解曲面の平均曲率が 0，左側で 1 ということであるから，この境界線を横切って 2 階微分値は不連続になる．境界線が滑らかならば，この線に沿って u は滑らかに変化し，もちろん領域 G 内の境界線を除いたところでは

u は滑らかに変化する．
楕円型微分方程式 (11.15) の解 u は，領域 G 内で f, aik , bi , c ∈ C∞ であれば，境界 Γ が折れ曲がり，あるいは境界値が不連続でも，双曲型微分方程式とは異なり，領域 G の内部では解析的である．数学ではこのようにいわれるが，数値解析ではこのような境界の近傍点には，境界の不連続性を考慮した差分式を用いることが望ましい．下図に示すように境界 Γ が点 P で折れ曲がりここで領域 G が凹になっている場合には，式 (11.15) のすべての係数と境界値が滑らかでも，解は折れ曲り点で連続微分可能にはならない．点 Pの近傍で関数 u(r, θ) はu(r, θ) = uP +a0 (θ)rα +a1 (θ)r2α +· · ·
π
α=
2π−φ

(1/2 ≤ α < 1)のような級数に展開され，半径 r に関し α 次の Hölder 連続になる．数値解析においてもこのことを考慮し，折れ曲り点とその近傍点に対しては特別の差分式を用意することが望ましい5 ．

楕円型微分方程式 (11.15) に関しては，領域 G 内で f, aik , bi , c ∈ C∞ で，その滑らかな境界の部分 Γ0 上に解析的に滑らかな境界値が与えられるときに，G+Γ で解析的な解 u は，Γ0 を越えてその外部領域へ解析接続することができる．すなわち G の近傍へ u の値を一意的に外挿することができるが，この初期値問題
はいわゆる不適切な問題 (improperly-posed problem) で，解は通常波打つことになり，平滑化 (smoothing)の操作が必要でる．またたとえ平滑化を行ったとしても大幅延長は困難である．

5 今 3 項目まで取ることにすれば u(r) = u + a r α + a r 2α ，ポテンシャル方程式の差分解法に必要な，折れ曲り点の隣接格
001
子点の 2 階の差分式は次のように導かれる．u00 (r) = a0 α(α − 1)r α−2 + a1 2α(2α − 1)r2α−2 ，未定係数 a0 , a1 は格子点条件α2αui = u0 + a0 (i∆r) + a1 (i∆r) , (i = 1, 2) から決定され，次の差分式が得られる．
i¯˘¯
1 1 h α˘
α2 (α−1)2α − 2(2α − 1) (u1 − u0 ) + α − (α−1) + 2(2α − 1) (u2 − u0 )
u00
1 =
∆r2 ∆

ただし ∆ = 2α (2α − 1)．

11

回転流面上の翼列流れ

11.2

回転流面上の２次元翼列流れの FEM 解析

軸流型，斜流型，遠心型ターボ機械の羽根車の流れを，１つの回転流面 (revolutional flow surface) を羽根車の流路内に取り，この回転流面上の２次元流れとして解析する方法について述べる．このような翼列流れの計算法は，コンピュータの未発達の時代に開発されたもので時代遅れの感は否めないが，２次元流れの有限要素法 (FEM, finite element methods) 解析の１つの例として見ていただければと思う．

11.2.1

流れ関数方程式の境界値問題

始めに軸流型の場合を考えれば，回転流面は半径 r の円筒面になり，その流れは z/r = ξ, θ = η と置くことによって，ξη 平面上の直線軸翼列を通る流れに写像して解析することができる．zθr は円柱座標である．この流れの方程式は次の Laplace 方程式になる．

∂2Ψ ∂2Ψ
+=0∂ξ 2 ∂η 2(11.16)

ただし Ψ は羽根車に固定された座標系から見た流れすなわち相対流れの流れ関数 (stream function) で，流速の z, θ 方向成分は次式で与えられる．

1 ∂Ψ
,
r ∂θ

wz =

wu = −∂Ψ∂z(11.17)次に一般的な斜流型の場合を考える．回転流面上に子午線方向と周方向の座標 mθ を取り，座標変換
dm/r = dξ, dθ = dη を行い，mθ 面上の翼列を 軸流型の場合と同様に ξη 面上の直線軸翼列に等角写像すれば，写像面上の流れの方程式は次の Poisson 方程式になる．

∂2Ψ ∂2Ψ
∂
∂Ψ
∂
∂Ψ
+ 2 =
(ln bρ)+ (ln ρ)+2ωr2 bρ cos δ
2
∂ξ
∂η
∂ξ
∂ξ ∂η
∂η

(11.18)

またこの場合の流れ関数は次式で定義される．

wm =

1 ∂Ψ
,
rbρ ∂θ

wu = −1 ∂Ψbρ ∂m(11.19)ただし b は回転流面に付加される流路幅方向の厚み，δ は回転流面の内向き法線と回転軸のなす角，ω は羽根車の回転の角速度，ρ は流体の密度である．なお同軸円筒壁を持つ軸流型の場合には b = 1, cos δ = 0，また非圧縮性流れの場合には ρ = 1 である．
b(ρwu +∆θ×∂(ρwu )/∂θ)∆m(bρwm r+∆m
×∂(bρwm r)/∂m)∆θ−bρwm r∆θ−bρwu ∆m

P
PP ¡
##
@
P¡
@#∆m
P
P
P
#
#
#
P
PPr∆θ P
#
PPP
#
@
@
¡
¡
(a)

−(wm +∆θ×∂wm/∂θ)∆m(wu r+∆m×∂(wu r)/∂m)∆θ

P
P
¡
@
#P
@## P
J¡
P
ω
−2ω
J
#
#
#
P
#H
PP J #
PP# H wm ∆mJ^#¡¡−wu r∆θ−2ω cos δ(b)

式 (11.18) は次のように導出される．(a) 図に示すように回転流面上に要素 ∆m×r∆θ を取り，この要素に対し相対流れの発散を取れば次の連続の式が得られる．

∂∂(rbρwm ) +(bρwu ) = 0
∂m
∂θ

12

非粘性流れの FEM 解析

式 (11.19) はこの式を満足するように定義された流れ関数 Ψ の式である．また (b) 図に示すように上記の要素を一巡する循環を取れば，相対流れ場には渦度 −2 ω が存在し，この要素に対する法線成分は −2ω cos δであるから，次式が得られる．

∂∂wm(rwu ) −= −2ωr cos δ∂m∂θ
この式に式 (11.19) を代入すれば∂ ³ r ∂Ψ ´
∂ ³ 1 ∂Ψ ´+
= 2ωr cos δ∂m bρ ∂m∂θ rbρ ∂θ

(11.20)

更に座標変換 r∂/∂m = ∂/∂ξ, ∂/∂θ = ∂/∂η を行えば式 (11.18) が導かれる6 ．
式 (11.18) の境界値問題では，回転流面の形状 r = r(z)，回転流面上の羽根の形状 θ = θ(m)，羽根車の羽根枚数 N ，回転の角速度 ω ，上流境界の子午線速度 wmU などが与えられる．ここでは遠心ポンプの羽根車を想定し，回転流面に流路高さ b を流路断面積一定すなわち子午線速度の流路断面にわたる平均値が一定になるように付加することにする．br = const. ≡ 1．このとき式 (11.18) と (11.19) は次のようになる．

Ψξξ + Ψηη = −(ln r)ξ Ψξ + 2ωr cos δ(11.21)

wm = Ψη ,(11.22)

wu = −Ψξ

なお式 (11.21) の右辺は次のように書換えることもできる．

rhs = −(ln r)ξ Ψξ + 2ωr cos δ = (wu +2ωr) cos δ(11.23)

ケーシングに固定された絶対座標系の絶対速度 v と羽根車に固定された相対座標系の相対速度 w の間には，

ω ×rr のような関係があり，これを図示すれば次の速度線図と呼ばれるものになる．この速度線図
v = w +ωから次の関係が得られる．
w2 = v 2 +(ωr)2 −2ωrvu−1

α = tan(11.24)
6
v >
 ¡
¡ vu


(vu /vm )β = tan−1 (wu /wm )

ω ×rr


@ vm = w m
@
@
@
w@
R
@

wu

流れ関数の方程式 (11.21) の境界値問題の境界条件は次のように与えられる．翼型表面には Dirichlet 条
件，上流と下流の境界には Neumann 条件，翼列の上流と下流領域の周期境界には周期境界条件が与えられる．
下方の翼面：

Ψ = 0,

上流と下流境界：
周期境界：

上方の翼面：

Ψ = q,

Ψξ = −wu ,Ψ (η0 +t) = Ψ (η0 )+q,

(11.25a)(11.25b)Ψη (η0 +t) = Ψη (η0 )6 u(x) に対し d(ln u)/dx = (1/u)du/dx，また ud(1/u)/dx = d{ln(1/u)}/dx = −d(ln u)/dx．

(11.25c)

13

回転流面上の翼列流れ

ただし q = wm t は１つの翼間流路の流量，t = 2π/N は ξη 面上の直線軸翼列のピッチである．また上流境
界の wu は与えられた vuU から直接求められるが，下流境界の wu は，Kutta の条件すなわち翼後縁から流れが滑らかに流下するという条件が満足されるように予測子修正子法によって求められる．
この境界値問題は通常 差分法または有限要素法 (FEM) で解かれる．ここではこの境界値問題を三角形 1
次要素を用い Galerkin FEM で解く．FEM については 9.1 節に一通り述べたので，ここには必要最低限を納得のいくように説明する．FEM で得られた Ψ の節点値の連立 1 次方程式は緩和法で解くことにする．

11.2.2

要素分割と補間関数

有限要素法では，まず閉じた計算領域を要素 (elements) に分割する (要素分割，mesh discretization)．下図に三角形 1 次要素と四辺形双 1 次要素を示す．• で示す点は節点 (nodal points) と呼ばれ，ここに座標と変数が定義される．メッシュは，要素が正三角形や長方形に近く，その大きさが徐々に変化し，精度の必要なところが細かくなっているものが望ましい．要素には要素番号，節点には節点番号，各要素の節点には要
素節点番号が振られる．配列 M を宣言し，あらかじめ各要素の節点番号を入れておく．
M( [要素節点番号], [要素番号] ) = [節点番号]
•3

1

¡@
¡
@
¡
@
Ã• 2
¡ ÃÃÃÃÃÃ¡•Ã
三角形 1 次要素Ã•
3 •ÃÃÃ
B

1

¯
•¯

¯
¯

4

B

B

B

B•

2

四辺形双 1 次要素

要素内の任意の点における変数の値は節点値から補間によって求められる．補間関数 (interpolation function)
は，最も簡単な三角形 1 次要素の場合には x, y の 1 次式で与えられる．すなわちx) = c1 +c2 x+c3 yu∗ (x(11.26)

ただし ci , (i = 1, 2, . . . ) は未定係数で，その数は一般に要素の節点数に等しくなければならない．また式(11.26) は節点上でも成り立つから，
ui = c1 +c2 xi +c3 yi(i = 1, 2, 3)(11.27)のような節点条件が得らる．これは ci の連立 1 次方程式で，これを解けば未定係数 ci したがって補間関数を決定することができる．しかしながら通常は補間関数は次の方法で求められる．その方が一般に容易に求められるからである．

u∗ は式 (11.26) により c1 , c2 , c3 の 1 次式，c1 , c2 , c3 は式 (11.27) により u1 , u2 , u3 の 1 次式，したがって u∗ は u1 , u2 , u3 の 1 次式である．補間関数 u∗ は一般に変数 u の節点値の 1 次式で次のように表すことができる．

x) = φ (xx) · uu∗ (x(11.28)x) は形状関数 (shape function，基底関数 base function) と呼ばれるもので，節点と補間しようと
ただし φ (xする点の座標のみの関数である．また u は変数 u の節点値のベクトルである．これらは三角形 1 次要素の14

非粘性流れの FEM 解析

場合には次のようになる．

φ = (φ1 φ2 φ3 ),

 
u1
 
u = u2 
u3(11.29)式 (11.28) に式 (11.26) と (11.27) を代入すれば

c1 +c2 x+c3 y = (c1 +c2 x1 +c3 y1 )φ1 +(c1 +c2 x2 +c3 y2 )φ2 +(c1 +c2 x3 +c3 y3 )φ3
この式は c1 , c2 , c3 の値によらず一般に成立すべきものであるから，c1 , c2 , c3 にかかる両辺の係数を等置す
れば次式が得られる．

1 = φ1 +φ2 +φ3x = x1 φ1 +x2 φ2 +x3 φ3y = y1 φ1 +y2 φ2 +y3 φ3
これは φi の連立 1 次方程式で，Cramer 法で解けば次の三角形 1 次要素の形状関数の式が得られる．
¯ ¯
¯ ¯
¯
¯
¯ 1 x y ¯ ¯ 1 x1 y1 ¯ ¯ 1 x1 y1 ¯
¯
¯
¯
¯
¯
¯
1 ¯
¯ ¯
¯ ¯
¯
φ = (φ1 φ2 φ3 ) =(11.30)¯ 1 x2 y2 ¯ ¯ 1 x y ¯ ¯ 1 x2 y2 ¯
¯ ¯
¯ ¯
¯
∆ ¯
¯ 1 x3 y3 ¯ ¯ 1 x3 y3 ¯ ¯ 1 x y ¯ただし¯¯ 1 x1¯¯∆ = ¯ 1 x2¯
¯ 1 x3
この行列式は
¯
¯x −x
¯ 2
1
=¯
¯x3 −x1

¯
y1 ¯¯
¯
y2 ¯
¯
y3 ¯

(11.31)¯¯
y2 −y1 ¯¯ ¯¯x2 −xx1 )×(xx3 −xx 1 )¯¯ = (x¯y3 −y1のように書換えることができるから，三角形要素面積の 2 倍であることが分かる．式 (11.30) の右辺の行列x) は，x 1 , x 2 , x 3 を頂点とする三角形の面積に対する式も同様に解釈することができ，これより例えば φ1 (xx) を示せx , x 2 , x 3 を頂点とする三角形の面積の割合であることが分かる．xy 面に垂直に φ 軸を取り φ1 (xば，節点 x 1 で高さ 1，節点 x 2 , x 3 で高さ 0 の平面になる．FEM では関数 u の微分積分は補間関数を導入
し各要素ごとにその形状関数を微分積分することによって計算される．例えば
Z
Z
∂φi
∂u∗
∗
=
ui ,
u dΩ =
φi dΩ ui
∂xj
∂xj
Ωe
Ωe
R
ただし Ω e dΩ は 1 つの要素にわたる積分である．三角形 1 次要素の場合には，φ は x, y の 1 次式であるからその 1 階微分は定数になり7 ，
Ã
Ã !
φ ,x
1 y2 −y3 y3 −y1=∆ x3 −x2 x1 −x3
φ ,yy1 −y2x2 −x1!
=1
Φ∆

7 添字 ‘, x’ は x に関する偏微分を意味する．φ,x = ∂φ/∂x．

(11.32)

15

回転流面上の翼列流れ

2 階微分は 0 になる．また φ の積分は結果だけ示せば次のようになる8 ．
Z
l! m! n!
(φ1 )l (φ2 )m (φ3 )n dΩ =∆(l+m+n)!
Ωe(11.33)四辺形双 1 次要素の補間関数も三角形 1 次要素と同様の方法で簡単に求めることができるが，ここでは
x 面上の四辺形要素を ξ 面上の正方形要素 (−1 ≤ ξ, η ≤ 1) に写像し，写像面上で補間する方法について述べる．この正方形要素の補間関数は，双 1 次すなわち ξ または η に関して 1 次で，また形状関数の性質上
例えば φ1 (ξξ i ) は i = 1 で 1，i = 2, 3, 4 で 0 でなければならないから次式で与えられる．

u∗ (ξξ ) = φ (ξξ ) · u = φ1 u1 +φ2 u2 +φ3 u3 +φ4 u4 ,φk =1
(1+ξk ξ)(1+ηk η)
4(11.34)(k = 1, 2, 3, 4)

ただし ξk , ηk の値は次の通りである．

3•

k

1

2

3

4

ξk
ηk

−1
−1

1
−1

−1
1

1
1

•4
η6

ξ

k = 1••2形状関数の ξ, η に関する微分は次のようになる．

φk,ξ =
あるいは
Ã !
φ ,ξ

φ ,η

1
ξk (1+ηk η),
4

1
=
4

Ã

−(1−η)−(1−ξ)

φk,η =

1(1+ξk ξ)ηk4

1−η−(1+ξ)−(1+η)
1−ξ

1+η
1+ξ

!

また形状関数の x, y に関する微分は
Ã ! Ã
!Ã !
Ã !
φ ,ξ
x,ξ y,ξ
φ ,x
φ ,x
=
=J
φ ,η
x,η y,η
φ ,yφ ,yを解くことによって求められる．J はヤコビ行列でアイソパラメトリック要素9 では x,ξ = φ ,ξ (ξξ )· x , · · · であるから次のようになる．

Ã !
´
φ ,ξ ³
J=
x y ,
φ ,η


x1
³
´ 
x2
x y =
x
 3
x4


y1

y2 

y3 

y4

なおここでは x は座標 x の節点値のベクトルである．また積分は例えばZ

Z 1Z 1x)φl (xx)dΩ =φm (xΩe

φm (ξξ )φl (ξξ )J(ξξ )dξdη =
−1 −1

4
X

J(ξξ k )φm (ξξ k )φl (ξξ k )(11.35)k=1となる．J = |J| はヤコビアン，積分値は通常 Gauss 積分 (6.1.3 項に述べた Gauss-Legendre 求積法) で計√算され，ξk , ηk はその積分点 (Gauss 点) の座標で上式では ±1/ 3 である．
8 9.3 節 有限要素法参照．

9 アイソパラメトリック要素 (iso-parametric element) は，関数 u と座標 x, y に対し同じ節点と形状関数を用いるものである．

16

非粘性流れの FEM 解析

11.2.3

ガレルキン法によるポアソン方程式の解法

FEM は変分 FEM，重み付き残差法，最小２乗 FEM などに分けられ，また重み付き残差法 (methods ofweighted residuals) は重み関数 (weighted functions) の選び方によっていくつかの方法に分けられる．ガレ
ルキン法 (Galerkin method) は重み関数に形状関数を用いるもので最も広く用いられている．Poisson 方程式の混合境界値問題は一般に次のように書くことができる．

∇2 u ≡ u,xx +u,yy = f (x, y)x∈Ω(11.36a)u = g(s)x ∈ ∂Ω1(11.36b)u,n = h(s)x ∈ ∂Ω2(11.36c)

ただし ∂Ω は領域 Ω の境界で ∂Ω = ∂Ω1 +∂Ω2 ，∂Ω1 は Dirichlet 条件すなわち関数 u の値の与えられる境界，∂Ω2 は Neumann 条件すなわち関数の法線微分 n ·∇u = u,n の値の与えられる境界である．なお，n は境界における外向き単位法線ベクトルである10 ．

Galerkin 法では，式 (11.36) を直接解くのではなく，この式に形状関数 φk を掛け計算領域に亙って積分した次式が解かれる．
ZZ(u,xx +u,yy −f )φk dΩ −Ω(u,n −h)φk ds = 0(k = 1, 2, . . . , K)

(11.37)∂Ω2

形状関数 φk は，xyφ 空間において k 番めの節点で高さ 1，この節点を共有する要素の集合の外側で高さ 0の角錐状のものになり，通常 三角形 1 次要素では六角錐になり，また四辺形双 1 次要素では稜線の大きく崩落したピラミッド状のものになる．節点が境界上にある場合にはこれらの錘体は中央で切断される．上式の積分は，φk 6= 0 の範囲，すなわち節点 k を共有する要素の集合に対して行われる．
Z
XZ
dΩ =
dΩ
Ω

n

Ωe

ただし n は要素番号で，P

n は形状関数の値が 0 でない要素の和を意味する．

ところで変分問題では，汎関数の第１変分を取るときに，Euler 方程式と共に 関数 u の境界値の与えられない境界の部分で自然境界条件 (natural boundary condition) が得られる．極値関数は Euler 方程式を
自然境界条件を考慮して解くことによって求められる．この汎関数の変分問題と微分方程式の境界値問題の類似性を考えれば，式 (11.37) のような Galerkin 法の式は第 1 変分をとるときに得られる式，微分方程式は Euler 方程式，Neumann 条件は自然境界条件，境界値問題の解は極値関数に対応する．FEM では，

Neumann 条件はこのような意味合いにおいて自然境界条件と呼ばれ，また u,n = 0 の場合に ‘自然に’ 扱われたと言うのは，何もしないのにひとりでに処理されたとも解釈できるが，自然境界条件として処理されたという意味であろう．微分方程式は，差分法の定式化では直接的に差分方程式で近似される．これに対し
FEM の定式化では，微分方程式は計算領域のある部分領域に亙って平均的に満足されるように処理される．
FEM のこのような定式化は，弱定式化 (weak formulation) と呼ばれる．
式 (11.37) は，部分積分することによって微分の階数を下げることができる．この式の第 1 項の被積分関数は(u,xx +u,yy )φk = φk∇2u = ∇·(φk∇u)−∇φk ·∇u = ∇·(φk∇u)−(φk,x u,x +φk,y u,y )10 数学では更に u(xx) が Ω +∂Ω で連続で，Ω で正則という条件が課される．第 1 の条件は領域内部の解が境界に繋がっているということ，第 2 の条件は領域内部で上記問題では 2 回連続微分可能 (twice continuously differentiable)，すなわち関数自身とその1 階と 2 階微分が Lipschitz 連続であるということである．

17

回転流面上の翼列流れとなり，またこの式の右辺第 1 項の積分は Gauss の定理により，
Z
Z
Z∇·(φk∇u)dΩ =n ·φk∇u ds =φk u,n ds
Ω

∂Ω

∂Ω2

となるから，これらの関係を式 (11.37) に用いれば次式が得られる．
ZZ(φk,x u,x +φk,y u,y +φk f )dΩ −φk h ds = 0Ω∂Ω2

なお h = 0 の場合には Neumann 条件は何もしなくてもひとりでに処理されることになり，これが FEM の利点の一つになっている．上式の u, f, h に式 (11.28) の補間関数を用いれば次式が得られる．

X

e
Kkj
uj = −

n

X

e
Fkj
fj +

n

X0

e
Hkj
hj

(k = 1, 2, . . . , K)

(11.38)

n

ただし

Z
e
Kkj
=

また

Z
(φk,x φj,x +φk,y φj,y )dΩ,

Ωe

e
Fkj
=

Z
φk φj dΩ,

Ωe

e
Hkj
=

φk φj ds

(11.39)∂Ω2e

P0

n は k 番め節点の隣接要素の，Neumann 条件の課されている辺の和である．

最終的に解かれる uj の連立 1 次方程式は式 (11.38) に Dirichlet 条件 (11.36b) を考慮したものになる．

Akj uj = bk(k = 1, 2, . . . , K)

A u = b)(A

(11.40)

ただし

(P
Akj =

e
n Kkj

δkj
 P
e

 − n Fkj fj
bk =
gk

P0 e
 P e
hj
− n Fkj fj + n Hkj

x k ∈ Ω +∂Ω2x k ∈ ∂Ω1

(11.41)xk ∈ Ω
x k ∈ ∂Ω1
x k ∈ ∂Ω2

(11.42)

δkj は Kronecker のデルタ関数，また x k ∈ Ω +∂Ω2 は k 番めの節点が領域内部の点か Neumann 条件の課される境界上の点であることを意味する．式 (11.40) の k 番めの式は k 番めの節点の形状関数を重み関数に用いた k 番めの節点の式で，この式全体は節点数 K に等しい未知数 uj を持つ K 個の式からなる連立 1次方程式である．式 (11.40) の係数行列 A を剛性マトリックス (stiffness matrix)，右辺の b を外力項 (forceterm) と言うことがあるが，これは FEM がはじめ構造解析の分野で発展し，A が剛性，b が外力を表していたからである．上式の A は対称でないが，式 (11.41)，(11.42) で
δkj −→Xe
Kkj+ λδkj ,n

gk −→ λgkのように置き λ を十分大きく取れば対称にすることができる．

18

非粘性流れの FEM 解析

11.2.4緩和法

連立 1 次方程式 (11.40) は直接法または緩和法で解かれる．ここには緩和法について述べる．翼列の上流
または下流領域では，式 (11.22) の第 2 式を周方向に平均すれば

Ψ ξ = −wu = ωr−Γ
2πr(11.43)となる．この式は計算点が ξ 一定線上に並びかつ η 方向に等間隔の場合には次のようになる．

¡
Γ ¢
∆η
q
(11.44)
Ψ i+1 = Ψ i + ωr−
∆ξ +
i+1/2
2πr
t
¢
¡ Pjf
ただし Ψ i =j=1 Ψij /jf , Γ = ΓU または ΓD , ∆ξ = ξi+1 −ξi , ∆η = ηi+1,0 −ηi0 である．翼列の上流と下流領域では，Ψ の予測値は式 (11.44) を基に与えられる．またこれらの領域の緩和計算では，予め式(11.44) を満足するように列ごとに一様な修正値を加えてから通常の緩和計算を行う．
点過緩和法では，各点ごとに式 (11.40) から残差
(n)

rk

(n)

= Akj uj

(n)

− bk

を計算し，この残差を打ち消すように次式によって修正計算が行われる．
(n+1)

uk

(n)

(n)

= uk −αrk /Akk

ただし α は緩和係数で 1.5 程度に取られる．下記のプログラムでは列ごとに一括して緩和計算を行う線過緩和法 (line over-relaxation method) が用いられる．上流と下流の境界では周方向に一様な流れを仮定しているので，ここでは通常の緩和計算は行わない．なおこれら境界の近傍に大きい残差が残るようであれば計算領域を広げることを検討しなければならない．

11.2.5

FEM における精度の良い数値微分

Laplace 方程式や Poisson 方程式の，三角形 1 次要素や四辺形双 1 次要素を用いる FEM 解析の精度は，これらの要素が正三角形または長方形またはこれらに近い場合には，2 次精度の 5 点差分式を用いる差分法に比べ良いことが知られている．しかしながらその微分値は，差分法では通常 2 次精度であるのに対し，

FEM では単に各要素内で補間関数を微分したのでは，2 次精度を得ることは難しい．以下には，これらの要素を用い FEM 解析で得られた結果から微分値を精度良く求める１つの方法について述べる．
•7

•3

1
•

\

6
•

•
1\

¿\
¿\
¿
¿
\
\
\• 5
¿
¿
◦
2•
\
¿
¿
¿0\
\
\
\•¿
\•¿
3

¿
¿

4

(a) 三角形 1 次要素の通常点

¿\

\
◦¿
0\
\

\
\•
¿5
¿

\• 6
\
¿\
¿
¿
\
\
\• 5
¿
¿
◦
1•
\
\
¿
¿0
¿
\
\•¿
\•¿
2\
4
¿
\
\•¿
8 •¿

•¿6

•
2

\•¿
4\

(b) 同境界点

図 11.1: 近傍節点の選び方¿\¿

3(c) 四辺形双 1 次要素

19

回転流面上の翼列流れのプログラムx0 ) = φ0 を満足する点 x 0 まわりの 2 次多項式は一般に次のように書くこと一つの節点 x 0 の節点条件 φ(xができる．

x ) = φ0 +
φ(x

L
X

cl ψl ( x )(11.45)l=1

ただし ２次元では L = 5, ψ1 = x, ψ2 = y, ψ3 = x2 , ψ4 = y2 , ψ5 = 2 xy である．x = x−x0 , y = y−y0 は
節点 x 0 に原点を持つ局所座標．式 (11.45) の L 個の未定係数 cl は，点 x 0 近傍の K 個の節点の条件

φk −φ0 =LX

cl ψlk(k = 1, . . . , K)

(11.46)l=1

が平均的に満足されるように，最小２乗法によって決定される．ただし K > L, φk = φ(xxk ), ψlk = ψl (xxk )である．すなわち cl の値は次の連立１次方程式を解くことによって求められる．
LX

anl cl = fn

(n = 1, 2, . . . , L)(11.47)

Wk 2 ψnk (φk −φ0 )(11.48)

l=1

ただし

anl =

K
X

Wk 2 ψnk ψlk ,

k=1

fn =

K
X
k=1

この式の Wk は重み関数である．微分値は次式で与えられる．

∂φ/∂x = c1 ,∂φ/∂y = c2

(11.49)

この方法で求められる微分値の精度は，重み関数と近傍点の選び方による．通常の最小２乗法では点 x 0
から遠い節点条件がむしろ重視される傾向にあるが，ここでは点 x 0 に近い節点条件を重視するために，重み関数を Wk 2 = (x2 +y2 )−3 のように置くことにする．また近傍節点は，図 11.1 に示すように，未定係数の数 L よりも多い点 x 0 を共有する要素のすべての節点を取り，また境界では領域の外にはみ出した点を内部に補う．2 次多項式近似を実現するには，このように近傍節点を点 x 0 のまわりのすべての方向に均等に，かつ 3 点ずつ並ぶように取ることが必要である．この方法によれば微分値は，適切に分割されたメッシュに対して，2 次相当の精度で求めることができる．

11.3

回転流面上の２次元翼列流れのプログラム

ここでは遠心ポンプの羽根車を通る流れを，回転流面上の翼列を通る定常２次元流れとして解析する．すなわち上記の流れ関数の Poisson 方程式の境界値問題を，Galerkin FEM で有限要素方程式 (連立 1 次方程式) に置き換え，これを緩和法で解く．そのプログラムは下記のようになる．メインプログラムでは，必要なデータを入力し，メッシュ分割，流れ関数の境界値問題の解法，流速と静圧の計算のサブルーチンを引用し，結果を出力する．i, j は回転流面上の計算点の子午線方向と周方向の番号で，i = 0, i1 , i2 , if は上流境界，羽根の前縁，後縁，下流境界，j = 0, jf は羽根の圧力面または負圧面とそれに繋がる周期境界である．また m は節点番号，n は要素番号．回転流面の形状は子午線上のいくつかの基準点の座標 zs , rs で与えられる．ls は基準点の子午線方向の番号 i を記録する配列．また羽根の形状は θi で与えられる．このほかに，omg：羽根車の回転角速度 ω ，nob：羽根枚数 N ，wmu, vuu, p0u：上流境界における子午線速度

20

非粘性流れの FEM 解析

wmU ，絶対流れの周方向速度 vuU ，岐点圧力 p0U = 1 の値が与えられる．これらは計算を始めるに際して必要十分なものである．上流側の循環 Gmu = ΓU は予め与えられた予旋回 (pre-rotation) の周速度 vuU から直接的に求められるが，下流側の循環 Gmd = ΓD は Kutta の条件を満足するように予測子修正子法によって求めなければならない．このために全体の計算も反復する必要がある．最後に節点の座標 z, θ, r，流れ関数 Ψ ，相対速度 (wm , wu )，静圧 p の計算結果をプリントし，また計算メッシュと静圧分布＋流線のグラフを描く．以下に示すプログラムは FORTRAN 95/90 Free source form で書かれている．
PROGRAM MAIN!*******************************************************************************************! Problem: 2D Incompressible Inviscid Flow through a Cascade on a Revolutional Flow Surface! Numerical Methods: Galerkin FEM Using Triangular Linear Elements, Relaxation Method!*******************************************************************************************DIMENSION z(0:30),th(0:30,2),r(0:30),xm(0:30),xi(2,217),zs(9),rs(9),ls(9),Psi(217), &res(217),w(2,217),nm(3,360),p(217),f(217)
COMMON /COND/omg,pi,nob,wmu,vuu,Gmu,p0u,tbt /IJP/i1,i2,ip,jp,na &/INDX/Gmd,DGmd,resmax,resmin
CHARACTER*8 z1,c(2)! ***** Input Geometries of Revolutional Flow Surface and Blade, and conditions
DATA i1,i2,if,jf,lf/10,26,30,6,9/
ip=if+1; jp=jf+1; mf=ip*jp; nf=2*if*jf
DATA zs/ 0., 13.0, 37.5, 57.0, 75.0, 90.3, 98.3, 2*99.0/
DATA rs/2*75.0, 76.0, 81.0, 89.5, 107.0, 135.0, 159.0, 185.0/
DATA ls/0,2,6,10,14,18,22,26,30/
DATA th/10*0., 2.200, 2.084, 1.962, 1.817, 1.662, 1.517, 1.364, 1.207, 1.050, .894, &.736, .592, .465, .350, .242, .126, 5*0.,&
10*0., 2.200, 1.993, 1.829, 1.680, 1.528, 1.380, 1.228, 1.072, .921, .763, &.607, .462, .333, .227, .145, .072, 5*0./
DATA omg,pi,nob,wmu,vuu,p0u/1., 3.14159, 6, .23, .04, 1./! omg:angular velocity, nob:number of blades, wmu:upstream meridian velocity! vuu:upstream circumferential absolute velocity, p0u:upstream stagnation pressure! ***** Discretize into elements
CALL SETXI(z,th,r,xm,if,jf,xi,mf,zs,rs,ls,lf)
na=0! ***** Compute stream function Psi and velocities W
100 na=na+1
CALL CALPSI(z,r,xm,if,jf,xi,Psi,res,w,mf,nm,nf)
CALL VELCTY(r,if,jf,xi,Psi,w,p,mf)! ***** Output computational results
OPEN(20,FILE=’OUTPUT.DAT’)WRITE(20,60)na,Gmd,DGmd,resmax,resmin! Gmd:downstream circulation, DGmd:inclement of Gmd, resmax and resmin:maximum or minimum! residual of finite element eqns
60 FORMAT(/1H ’na =’, I3, 5X ’Gmd =’, F8.4, 5X ’DGmd =’, F8.4, &5X ’resmax =’, F8.4, 5X ’resmin =’, F8.4)
IF(na<99)
GOTO 100
DO k=1,5
SELECT CASE(k)CASE(1); z1=’
Psi ’; isf=0; FORALL(m=1:mf)f(m)=Psi(m)CASE(2); z1=’ res*100’; isf=2; FORALL(m=1:mf)f(m)=res(m)CASE(3); z1=’
w_m ’; isf=0; FORALL(m=1:mf)f(m)=w(1,m)
CASE(4); z1=’
w_u ’; isf=0; FORALL(m=1:mf)f(m)=w(2,m)
CASE(5); z1=’pressure’; isf=0; FORALL(m=1:mf)f(m)=p(m)ENDSELECTCALL WRTF(f,z1,ifs)
ENDDO

回転流面上の翼列流れのプログラム
IF(na<100)
GOTO 100
CLOSE(20)! ***** Drawing mesh, pressure distribution and streamlines
CALL RESULTS(xi,Psi,p,if,jf,mf,nm,nf)
STOP
END PROGRAM MAIN
! ********** Mesh Discretization
SUBROUTINE SETXI(z,th,r,xm,if,jf,xi,mf,zs,rs,ls,lf)
DIMENSION z(0:if),th(0:if,2),r(0:if),xm(0:if),xi(2,mf),zs(lf),rs(lf),ls(lf), &xs(9),us(9),x1(0:30),u1(0:30),f(217)
COMMON /COND/omg,pi,nob,wmu,vuu,Gmu,p0u,tbt /IJP/i1,i2,ip,jp /INDX/Gmd,DGmd
CHARACTER*8 z1,c(2)pit=2.*pi/nob; Dpit=pit/FLOAT(jf); r2=rs(8)! z,theta,r:cylindrical coords, xm:meridian coord, xi=(xi,eta):mapping plane coords! ***** Generate revolutional flow surface by interpolation using exponetial spline
DO n=1,3
IF(n==1)THENFORALL(i=0:if)x1(i)=FLOAT(i); us(1)=0.
ENDIF
FORALL(l=1:lf)xs(l)=x1(ls(l))
SELECT CASE(n)CASE(1)
DO l=2,lf
dz=zs(l)-zs(l-1); dr=rs(l)-rs(l-1)
us(l)=us(l-1)+SQRT(dz*dz+dr*dr)/r2!normalization
ENDDO
CASE(2); FORALL(l=1:lf)us(l)=zs(l)/r2!normalization
CASE(3); FORALL(l=1:lf)us(l)=rs(l)/r2!normalization
ENDSELECT
dx=x1(2)-x1(0);IF(n==1)dx=2.
g1=(us(2)-us(1))/dx
dx=x1(30)-x1(26); IF(n==1)dx=4.
gf=(us(9)-us(8))/dx
CALL ESPL(xs,us,ls,lf,g1,gf,x1,u1,if)
SELECT CASE(n)CASE(1); FORALL(i=0:if)x1(i)=u1(i)CASE(2); FORALL(i=0:if) z(i)=u1(i)CASE(3); FORALL(i=0:if) r(i)=u1(i)ENDSELECT
ENDDO
xm(0)=0.
DO i=1,if
dz=z(i)-z(i-1); dr=r(i)-r(i-1)xm(i)=xm(i-1)+SQRT(dz*dz+dr*dr)!more exact meridian coords
ENDDO
xi(1,1)=0.
DO i=1,if; m=i*jp+1
dm=xm(i)-xm(i-1); rm=(r(i-1)+r(i))/2.
xi(1,m)=xi(1,m-jp)+dm/rm
ENDDO! ***** Determine blade outlet angle tan(beta_t) and circulation Gmd
i=i2; m=i*jp+1
Dxi=xi(1,m)-xi(1,m-jp)
bt1=ATAN((th(i,1)-th(i-1,1))/Dxi)
bt2=ATAN((th(i,2)-th(i-1,2))/Dxi)
tbt=TAN((bt1+bt2)/2.)!blade outlet angle

21

22

非粘性流れの FEM 解析

Gmd=2.*pi*r(i)*(wmu*tbt+omg*r(i))! ***** Genarate periodic boundaries
Gmu=2.*pi*r(0)*vuu
DO i=i1-1,0,-1; m=i*jp+1
rm=(r(i)+r(i+1))/2.; Dxi=xi(1,m+jp)-xi(1,m)
wu=Gmu/(2.*pi*rm)-omg*rm; Dth=Dxi*wu/wmu
th(i,1)=th(i+1,1)-Dth
ENDDO
DO i=i2+1,if; m=i*jp+1
rm=(r(i-1)+r(i))/2.; Dxi=xi(1,m)-xi(1,m-jp)
wu=Gmd/(2.*pi*rm)-omg*rm; Dth=Dxi*wu/wmu
th(i,1)=th(i-1,1)+Dth
ENDDO! ***** Discretize computational domain into elements
DO i=0,if; m0=i*jp+1
IF(i<=i1.OR.i>=i2)th(i,2)=th(i,1)
th(i,2)=th(i,2)+pit
DO j=0,jf; m=i*jp+j+1; t=j/FLOAT(jf)xi(1,m)=xi(1,m0)
xi(2,m)=(1.-t)*th(i,1)+t*th(i,2)
ENDDO
ENDDO
! ***** Print xm, z, r, theta, xi,eta
OPEN(20,FILE=’OUTPUT.DAT’)WRITE(20,60)xm,z,r,th
60 FORMAT(/1H ’ m=’, 31F8.3/1H ’ z=’, 31F8.3/1H/1H ’theta=’, 31F8.3/1H 6X
31F8.3/)
DO i=1,2
FORALL(m=1:mf)f(m)=xi(i,m)
z1=c(i); isf=0; CALL WRTF(f,z1,isf)
ENDDO
DATA c/’xi’,’ eta’/
END SUBROUTINE SETXI!downstr circulation!upstr circulation!upstream region bound!wmu=const!downstream region bound

!xi
!eta

’ r

=’, 31F8.3 &! *********** Calculation of Cascade Flow by FEM
SUBROUTINE CALPSI(z,r,xm,if,jf,xi,Psi,res,w,mf,nm,nf)
DIMENSION z(0:if),r(0:if),xm(0:if),xi(2,mf),Psi(mf),res(217),w(2,mf),nm(3,nf), &A(217,-8:8),B(217),f(217),cosdel(0:30),phi(2,3),c(6,3),rhs(6)COMMON /COND/omg,pi,nob,wmu,vuu,Gmu,p0u,tbt /IJP/i1,i2,ip,jp,na &/INDX/Gmd,DGmd,resmax,resmin
CHARACTER*8 z1
pit=2.*pi/nob; Dpit=pit/FLOAT(jf); flw=pit*wmu!flw:flow through a channel
i0=if-1; j0=jf-1; lf=jp+1; n1=jf
DATA al1,al2/.6, 1.5/
IF(na/=1)
GOTO 200! ***** Setup nodal point numbers of each elements
n=1; m=1
201 nm(1,n)=m; nm(2,n)=m+jp; nm(3,n)=m+jp+1
n=n+1; nm(1,n)=m; nm(2,n)=m+jp+1; nm(3,n)=m+1
n=n+1; m=m+1; IF(MOD(n,2*jf)==1)m=m+1
IF(n<=nf)
GOTO 201! ***** Predict values of Psi
DO i=0,if; m0=i*jp+1; DO j=0,jf; m=i*jp+j+1!Psi=0 on pressure side
Psi(m)=Psi(m0)+flw*j/FLOAT(jf)!Psi=flw on suction side
ENDDO; ENDDO! ***** Predict upstream and downstream velocities
FORALL(m=1:jp)

回転流面上の翼列流れのプログラム
w(1,m)=wmu; w(2,m)=vuu-omg*r(0)ENDFORALLFORALL(m=mf-jf:mf)
w(1,m)=wmu; w(2,m)=Gmd/(2.*pi*r(if))-omg*r(if)
ENDFORALL!at upstr bound!at downstr bound

GOTO 210! ***** correct Gmd and w_u at downstream bound200 m=i2*jp+1; m0=m+jp!at blade trailing edge
bwm=(w(1,m)+w(1,m+jf)+2.*w(1,m0))/4.
!b*=bar(*):averaged *bwu=(w(2,m)+w(2,m+jf)+2.*w(2,m0))/4.
btb=(tbt+w(2,m0)/w(1,m0))/2.
Dxi=xi(1,m0)-xi(1,m); Deta=xi(2,m)-xi(2,m0)
DGmd=al1*4.*pi*r(i2)*(bwu+bwm*(Deta+Psi(m0)/w(1,m0))/Dxi) !al1:damping factor
Gmd=Gmd+DGmd!corrected Gmd
FORALL(m=mf-jf:mf)w(2,m)=Gmd/(2.*pi*r(if))-omg*r(if)!at downstr bound
210 CONTINUE! ***** Compute rhs of Poisson eq of streamfunction Psi
DO i=1,if-1
dr=r(i+1)-r(i-1); dm=xm(i+1)-xm(i-1)
cosdel(i)=dr/dm!cosdel=cos(delta)
ENDDO
dr=r(1)-r(0);dm=xm(1)-xm(0);cosdel(0)=dr/dm
dr=r(if)-r(i0); dm=xm(if)-xm(i0); cosdel(if)=dr/dm
DO i=0,IF; DO j=0,jf; m=i*jp+j+1
IF(na==1)THEN
IF(i<=i1)THEN;
Gmg=Gmu
ELSEIF(i>=i2)THEN; Gmg=Gmd
ELSE; Gmg=(Gmu*(i2-i)+Gmd*(i-i1))/(i2-i1); ENDIF
w(2,m)=Gmg/(2.*pi*r(i))-omg*r(i)ENDIF
f(m)=(w(2,m)+2.*omg*r(i))*cosdel(i)!rhs of Poisson eqn on xi
ENDDO; ENDDO! ***** Set up coefficient matrix A and rhs B of linear equations
FORALL(m=1:mf,l=-lf:lf)A(m,l)=0.; FORALL(m=1:mf)B(m)=0.
DO n=1,nfm1=nm(1,n); m2=nm(2,n); m3=nm(3,n)
phi(1,1)=xi(2,m2)-xi(2,m3); phi(2,1)=xi(1,m3)-xi(1,m2)
phi(1,2)=xi(2,m3)-xi(2,m1); phi(2,2)=xi(1,m1)-xi(1,m3)
phi(1,3)=xi(2,m1)-xi(2,m2); phi(2,3)=xi(1,m2)-xi(1,m1)
AJ=phi(2,3)*phi(1,2)-phi(2,2)*phi(1,3)!JacobianB(m1)=B(m1)-(2.*f(m1)+f(m2)+f(m3))*AJ/24.
!rhs of finite element eqns
B(m2)=B(m2)-(2.*f(m2)+f(m3)+f(m1))*AJ/24.
B(m3)=B(m3)-(2.*f(m3)+f(m1)+f(m2))*AJ/24.
DO me=1,3; m=nm(me,n)
DO le=1,3; l=nm(le,n)-m
Ae=(phi(1,me)*phi(1,le)+phi(2,me)*phi(2,le))/2./AJ
A(m,l)=A(m,l)+Ae!coefs of finite element eqns
ENDDO; ENDDO
ENDDO! ***** Solve simultaneouse linear qns of Psi by relaxation method! ***** Compute cascade region by usual SLOR
DO ii=i1,i2; i=ii; IF(MOD(na,2)==0)i=i1+i2-ii!alternate sweep direction
DO j=1,j0; m=i*jp+j+1
c(J,1)=
A(m,-1)c(J,2)=al2*A(m,0)
c(J,3)=
A(m,1)
rhs(J)=B(m)-A(m,-lf)*Psi(m-lf)-A(m,-jp)*Psi(m-jp) &

23

24

非粘性流れの FEM 解析-(1.-al2)*A(m,0)*Psi(m)-A(m,jp)*Psi(m+jp)-A(m,lf)*Psi(m+lf)
ENDDO; rhs(j0)=rhs(j0)-A(m,1)*flw
n=j0
CALL GAUSS3(c,rhs,n1,n)
FORALL(j=1:j0)Psi(i*jp+j+1)=rhs(j)ENDDO! ***** Compute upstream region
Psim=0.; DO m=i1*jp+2,i1*jp+jp; Psim=Psim+Psi(m); ENDDO!add uniform correction
Psim=Psim/FLOAT(jf)DO i=i1-1,0,-1; m=i*jp+1
Dxi =xi(1,m+jp)-xi(1,m)
Deta=xi(2,m+jp)-xi(2,m); rm=(r(i)+r(i+1))/2.
Psim=Psim-(omg*rm-Gmu/(2.*pi*rm))*Dxi-flw*Deta/pit
Psia=0.; DO m=i*jp+2,i*jp+jp; Psia=Psia+Psi(m); ENDDO
Psia=Psia/FLOAT(jf); DPsi=Psia-Psim
FORALL(m=i*jp+1:i*jp+jp)Psi(m)=Psi(m)-DPsi
ENDDO
DO i=i1-1,1,-1!SLOR considered periodic cond
DO j=1,j0; m=i*jp+j+1
c(j,1)=
A(m,-1)c(j,2)=al2*A(m,0)
c(j,3)=
A(m,1)
rhs(j)=B(m)-A(m,-lf)*Psi(m-lf)-A(m,-jp)*Psi(m-jp) &-(1.-al2)*A(m,0)*Psi(m)-A(m,jp)*Psi(m+jp)-A(m,lf)*Psi(m+lf)
ENDDO
m=i*jp+2; rhs(1)=rhs(1)+c(1,1)*flw
m=i*jp+jp; c(jf,1)=A(m,-1); c(jf,2)=al2*(A(m,0)+A(m-jf,0)); c(jf,3)=A(m-jf,1)
rhs(jf)=B(m)+B(m-jf)-A(m,-lf)*Psi(m-lf)-(A(m,-jp)+A(m-jf,-jp))*Psi(m-jp) &-(1.-al2)*(A(m,0)+A(m-jf,0))*Psi(m)-(A(m,jp)+A(m-jf,jp))*Psi(m+jp)&-A(m-jf,lf)*(Psi(m+2)+flw)-c(jf,3)*flw
n=jf
CALL GAUSSP(c,rhs,n)
FORALL(j=1:jf)Psi(i*jp+j+1)=rhs(j)m=i*jp+1; Psi(m)=Psi(m+jf)-flw
ENDDO! ***** Compute downstream region
Psim=0.; DO m=i2*jp+2,i2*jp+jp; Psim=Psim+Psi(m); ENDDO!add uniform correction
Psim=Psim/FLOAT(jf)DO i=i2+1,if; m=i*jp+1
Dxi =xi(1,m)-xi(1,m-jp)
Deta=xi(2,m)-xi(2,m-jp); rm=(r(i-1)+r(i))/2.
Psim=Psim+(omg*rm-Gmd/(2.*pi*rm))*Dxi+flw*Deta/pit
Psia=0.; DO m=i*jp+2,i*jp+jp; Psia=Psia+Psi(m); ENDDO
Psia=Psia/FLOAT(jf); DPsi=Psia-Psim
FORALL(m=i*jp+1:i*jp+jp)Psi(m)=Psi(m)-DPsi
ENDDO
DO i=i2+1,if-1!SLOR considered periodic cond
DO j=1,j0; m=i*jp+j+1
c(j,1)=
A(m,-1)c(j,2)=al2*A(m,0)
c(j,3)=
A(m,1)
RHS(J)=B(m)-A(m,-lf)*Psi(m-lf)-A(m,-jp)*Psi(m-jp) &-(1.-al2)*A(m,0)*Psi(m)-A(m,jp)*Psi(m+jp)-A(m,lf)*Psi(m+lf)
ENDDO
m=i*jp+2; rhs(1)=rhs(1)+c(1,1)*flw
m=i*jp+jp; c(jf,1)=A(m,-1); c(jf,2)=al2*(A(m,0)+A(m-jf,0)); c(jf,3)=A(m-jf,1)
rhs(jf)=B(m)+B(m-jf)-A(m,-lf)*Psi(m-lf)-(A(m,-jp)+A(m-jf,-jp))*Psi(m-jp) &回転流面上の翼列流れのプログラム-(1.-al2)*(A(m,0)+A(m-jf,0))*Psi(m)-(A(m,jp)+A(m-jf,jp))*Psi(m+jp)&-A(m-jf,lf)*(Psi(m+2)+flw)-c(jf,3)*flw
n=jf
CALL GAUSSP(c,rhs,n)
FORALL(j=1:jf)Psi(i*jp+j+1)=rhs(j)m=i*jp+1; Psi(m)=Psi(m+jf)-flw
ENDDO! ***** Estimate residuals of finite element equations
DO i=1,if-1; DO j=1,jf-1; m=i*jp+j+1
res(m)=-B(m)+A(m,-lf)*Psi(m-lf)+A(m,-jp)*Psi(m-jp)+A(m,-1)*Psi(m-1) &+A(m,0)*Psi(m)+A(m,1)*Psi(m+1)+A(m,jp)*Psi(m+jp)+A(m,lf)*Psi(m+lf)
ENDDO; ENDDO
cycle_1: DO i=1,if-1; m=i*jp+j+1; IF(i>=i1.AND.i<=i2)CYCLE cycle_1!periodic boundj=0; m=i*jp+1; m1=(i+1)*jp
res(m)=-B(m)+A(m,-jp)*Psi(m-jp)+A(m,0)*Psi(m)+A(m,1)*Psi(m+1) &+A(m,jp)*Psi(m+jp)+A(m,lf)*Psi(m+lf) &-B(m1)+A(m1,-lf)*(Psi(m1-lf)-flw)+A(m1,-jp)*(Psi(m1-jp)-flw) &+A(m1,-1)*(Psi(m1-1)-flw)+A(m1,0)*(Psi(m1)-flw)+A(m1,jp)*(Psi(m1+jp)-flw)
ENDDO cycle_1
DO j=0,j0; m=j+1; mm1=m-1; IF(j==0)mm1=m!upstr boundres(m)=-B(m)+A(m,-1)*Psi(mm1)+A(m,0)*Psi(m)+A(m,1)*Psi(m+1) &!A(1,-1)=0+A(m,jp)*Psi(m+jp)+A(m,lf)*Psi(m+lf)-Dpit*w(2,1)
ENDDO
res(1)=res(1)-B(jp)+A(jp,-1)*(Psi(jp-1)-flw)+A(jp,0)*(Psi(jp)-flw) &+A(jp,jp)*(Psi(jp+jp)-flw)
DO j=1,jf; m=if*jp+j+1; mp1=m+1; IF(j==jf)mp1=m!downstr boundres(m)=-B(m)+A(m,-lf)*Psi(m-lf)+A(m,-jp)*Psi(m-jp)+A(m,-1)*Psi(m-1) &+A(m,0)*Psi(m)+A(m,1)*Psi(mp1)+Dpit*w(2,mf)!A(mf,1)=0
ENDDO
m=mf; m0=mf-jf
res(m)=res(m)-B(m0)+A(m0,-jp)*(Psi(m0-jp)+flw)+A(m0,0)*(Psi(m0)+flw) &+A(m0,1)*(Psi(m0+1)+flw)
resmax=MAXVAL(res,DIM=1)
resmin=MINVAL(res,DIM=1)
END SUBROUTINE CALPSI! ********** Compute velosities
SUBROUTINE VELCTY(r,if,jf,xi,Psi,w,p,mf)
DIMENSION r(0:if),xi(2,mf),Psi(mf),w(2,mf),p(mf),neb(6,217),bp(6),ps(5,6), &a(5,5),aw(5),iw(5)COMMON /COND/omg,pi,nob,wmu,vuu,Gmu,p0u,tbt /IJP/i1,i2,ip,jp,na
pit=2.*pi/nob; Dpit=pit/FLOAT(jf); flw=pit*wmu
rot=p0u-omg*Gmu/(2.*pi)!rothalpy! ***** Determine neighbouring nodal numbers of each nodes
CALL SETNEB(neb,if,jf,mf)! ***** Compute relative velocities w_m, w_u on revolutional flow surface
DO m=1,mf; m0=mi=(m-1)/jp; j=MOD(m-1,jp); IF(i==0.OR.i==if)
GOTO 301 !upstr & downstr bounds
DO k=1,6; mk=neb(k,m0)
bp(k)=Psi(mk)-Psi(m0)ps(1,k)=xi(1,mk)-xi(1,m0)
ps(2,k)=xi(2,mk)-xi(2,m0)
IF(i>=i1.AND.i<=i2)
GOTO 302 !cascade region
IF(j==0)THEN!periodic bound j=0
bp(2)=bp(2)-flw; ps(2,2)=ps(2,2)-pit
bp(3)=bp(3)-flw; ps(2,3)=ps(2,3)-pit; ENDIF
IF(j==jf)THEN!periodic bound jf

25

26

非粘性流れの FEM 解析

bp(5)=bp(5)+flw; ps(2,5)=ps(2,5)+pit
bp(6)=bp(6)+flw; ps(2,6)=ps(2,6)+pit; ENDIF
302 CONTINUE
ps(3,k)=ps(1,k)*ps(1,k)
ps(4,k)=ps(2,k)*ps(2,k)
ps(5,k)=2.*ps(1,k)*ps(2,k)
ENDDO
FORALL(n=1:5,l=1:5)a(n,l)=0.
DO n=1,5; DO l=1,5; DO k=1,6
wt=ps(3,k)+ps(4,k); wt=1./wt/wt/wt!weight functiona(n,l)=a(n,l)+wt*ps(n,k)*ps(l,k)
ENDDO; ENDDO; ENDDO
nf=5; e=1.E-5
CALL MATINV(a,aw,iw,nf,det,e)w(1,m)=0.; w(2,m)=0.
DO n=1,5; f=0.
DO k=1,6
wt=ps(3,k)+ps(4,k); wt=1./wt/wt/wt
f=f+wt*ps(n,k)*bp(k)ENDDO
w(1,m)=w(1,m)+a(2,n)*f!w_m
w(2,m)=w(2,m)-a(1,n)*f
!w_u
ENDDO
! ***** Compute static pressure p
301 CONTINUE
ww=w(1,m)*w(1,m)+w(2,m)*w(2,m)p(m)=rot+omg*omg*r(i)*r(i)/2.-ww/2.
!static pressure
ENDDO! ***** Take averaged values of P at leading and trailing edges
m=i1*jp+1; p(m)=.5*(p(m)+p(m+jf)); p(m+jf)=p(m)m=i2*jp+1; p(m)=.5*(p(m)+p(m+jf)); p(m+jf)=p(m)END SUBROUTINE VELCTY! ********** SUBSIDIARY SUBROUTINES! ********** Interpolation Using Exponential Spline
SUBROUTINE ESPL(xs,us,ls,lf,g1,gf,x,u,if)
DIMENSION xs(lf),us(lf),ls(lf),x(0:if),u(0:if),Dxs(9),a(9,3),b(9),c(9),p(9),d(9),e(9)FORALL(l=1:lf-1)
Dxs(l)=xs(l+1)-xs(l); p(l)=.1!Dxs=Delta x_s, _s:data points, p:tension of beam
ENDFORALL
n=0; 13 n=n+1
DO l=1,lf
FORALL(j=1:3)a(l,j)=0.; b(l)=0.
IF(l<lf)THEN
xp=Dxs(l)*p(l); p1=1./p(l)d(l)=-p1/xp+p1/TANH(xp)!de(l)= p1/xp-p1/SINH(xp)
!e
ENDIF
ENDDO
DO 10 l=1,lf; IF(l==1)
GOTO 11a(l,1)=e(l-1); a(l,2)=d(l-1)!coefs of lhs
b(l)=-(us(l)-us(l-1))/Dxs(l-1); IF(l==lf)
GOTO 10!rhs
11 a(l,2)=a(l,2)+d(l); a(l,3)=a(l,3)+e(l)!coefs of lhs
b(l)=b(l)+(us(l+1)-us(l))/Dxs(l)!rhs
c(l)=b(l)10 CONTINUE
b(1)=b(1)-g1; b(lf)=b(lf)+gf!g1,gf:gradients at both edges

回転流面上の翼列流れのプログラム
n1=9; n=lf
CALL GAUSS3(a,b,n1,n)ip=0
DO 12 l=2,lf-1
b0=b(l); c0=2.*c(l)/(Dxs(l-1)+Dxs(l))
IF(b0*c0>=0..AND..8*ABS(b0)<ABS(c0))
GOTO 12p(l-1)=p(l-1)+.1*n*n; p(l)=p(l)+.1*n*nip=ip+1
12 CONTINUE
IF(n<2.AND.ip>0)
GOTO 13l=1
DO i=0,if-1; IF(i==ls(l+1))l=l+1
al=(x(i)-xs(l))/Dxs(l); al1=1.-al
xp=Dxs(l)*p(l); p2=1./p(l)/p(l)u(i)=(us(l)-p2*b(l))*al1+(us(l+1)-p2*b(l+1))*al &+p2/SINH(xp)*(b(l)*SINH(al1*xp)+b(l+1)*SINH(al*xp))
ENDDO
u(if)=us(lf)END!compt u’’!spline & numer curvs!increase tension!al=alpha!compute interporated values! ********** Solution of Simultaneous Linear Eqns having Tri-Diagonal Reduced Array!
by Gaussian Elimination
SUBROUTINE GAUSS3(a,b,n1,n)
REAL a(n1,3),b(n1)DO k=1,n-1b(k)=b(k)/a(k,2); a(k,3)=a(k,3)/a(k,2)
b(k+1)=b(k+1)-a(k+1,1)*b(k); a(k+1,2)=a(k+1,2)-a(k+1,1)*a(k,3)ENDDOb(n)=b(n)/a(n,2)
DO k=n-1,1,-1b(k)=b(k)-a(k,3)*b(k+1)
ENDDO
END
! ********** Solution of Simultaneous Linear Eqns having Tri-Diagonal Reduced Array!
Modified by Periodic Bound Cond, by Gaussian Elimination
SUBROUTINE GAUSSP(a,b,n)
DIMENSION a(n,3),b(n)DO k=1,n-1b(k)=b(k)/a(k,2); a(k,3)=a(k,3)/a(k,2); a(k,1)=a(k,1)/a(k,2)
b(k+1)=b(k+1)-a(k+1,1)*b(k); a(k+1,2)=a(k+1,2)-a(k+1,1)*a(k,3)
a(k+1,1)=-a(k+1,1)*a(k,1)b(n)=b(n)-a(n,3)*b(k); a(n,2)=a(n,2)-a(n,3)*a(k,1); a(n,3)=-a(n,3)*a(k,3)ENDDOb(n)=b(n)/(a(n,1)+a(n,2)+a(n,3))
DO k=n-1,1,-1b(k)=b(k)-a(k,3)*b(k+1)-a(k,1)*b(n)
ENDDO
END
! ********** Input of Neighbouring Nodal Point Numbers
SUBROUTINE SETNEB(neb,if,jf,mf)
INTEGER neb(6,mf)
COMMON /IJP/i1,i2,ip,jp,na
DO k=1,6
SELECT CASE(k)CASE(1,2); ii=-1; CASE(3,6); ii=0; CASE DEFAULT; ii=1; ENDSELECT

27

28

非粘性流れの FEM 解析

SELECT CASE(k)CASE(2,3); ij=-1; CASE(1,4); ij=0; CASE DEFAULT; ij=1; ENDSELECT
DO m=1+jp,mf-jp
neb(k,m)=m+ii*jp+ij
ENDDO
ENDDO
DO i=i1,i2!on blade
m=i*jp+1; neb(2,m)=m+2*jp+2; neb(3,m)=m+2!press side
m=m+jf;neb(5,m)=m-2*jp-2; neb(6,m)=m-2!suction side
ENDDO
cycle_1: DO i=1,if-1; IF(i>=i1.AND.i<=i2)CYCLE cycle_1!periodic boundm=i*jp+1; neb(2,m)=m-2; neb(3,m)=m+jf-1!j=0
m=m+jf;neb(5,m)=m+2; neb(6,m)=m-jf+1!jf
ENDDO cycle_1
END! ********** Print of 2D Data on Revolutional Flow Surface
SUBROUTINE WRTF(f,z1,isf)
DIMENSION f(217),ff(0:30,0:6)
COMMON /IJP/i1,i2,ip,jp,na
CHARACTER*8 z1; CHARACTER*4 Form(7),Area(7)DATA form/’(1H ’,’2X I’,’2,4X’,’ 31(’,’’,’F8.4’,’)) ’/DATA area/’ -3P’,’ -2P’,’ -1P’,’’,’ 1P’,’ 2P’,’ 3P’/
OPEN(20,FILE=’OUTPUT.DAT’)WRITE(20,60)z1,naForm(5)=Area(isf+4)
FORALL(i=0:30,j=0:6)ff(i,j)=f(i*jp+j+1)
DO j=6,0,-1; WRITE(20,Form)j,(ff(i,j),i=0,30); ENDDO
WRITE(20,61)(i,i=0,30)
60 FORMAT(//1H 10X ’* * * ’, A8, ’ * * *’, 10X ’na =’, I4/)
61 FORMAT(1H 2X ’j/i=’, 31I8)
END! ********** Matrix Inversion by Gauss-Jordan Reduction
SUBROUTINE MATINV(a,aw,iw,n,det,e)
DIMENSION a(n,n),aw(n),iw(n)det=1.; e=AMAX1(e,1.E-5)
FORALL(k=1:n)iw(k)=k
DO k=1,nw=a(k,k); IF(ABS(w)>e)
GOTO 11l=k; 12 l=l+1; IF(l>n)STOP 5555
w=a(l,k); IF(ABS(w)<=e)
GOTO 12
FORALL(j=1:n)
aw(j)=a(k,j); a(k,j)=a(l,j); a(l,j)=aw(j)ENDFORALL
iw0=iw(k); iw(k)=iw(l); iw(l)=iw0
det=-det; 11 det=w*det
a(k,k)=1.
FORALL(j=1:n)a(k,j)=a(k,j)/w
cycle_1: DO i=1,n; IF(i==k)CYCLE cycle_1
w=a(i,k); a(i,k)=0.
FORALL(j=1:n)a(i,j)=a(i,j)-w*a(k,j)
ENDDO cycle_1
ENDDO
cycle_2: DO j=1,n; l=j; IF(iw(j)==j)CYCLE cycle_2
13 l=l+1; IF(iw(l)/=j)
GOTO 13
FORALL(i=1:n)回転流面上の翼列流れのプログラム

29

aw(i)=a(i,j); a(i,j)=a(i,l); a(i,l)=aw(i)ENDFORALL
iw(l)=iw(j)ENDDO cycle_2
END

サブルーチン CALXI： 計算領域をメッシュに分割する．回転流面の子午線はいくつかの基準点の座標 zs , rsでのみ与えられるので，始めに子午線の詳細な形状 zi , ri がこれらの座標を補間することによって求められる．その計算手順はまず m(i) 曲線上の点 mi を補間し，次に z(m) 曲線上の点 zi と r(m) 曲線上の点 ri を補間するというものである．これらの補間にはサブルーチン ESPL が用いられるがここでは指数スプラインをもちいて補間が行われる．g1,gf は指数スプラインの両端の勾配，x,u は補間点の座標 xi と関数値 ui の配列である．なお x1 はこの段階で用いられるおおよその子午線座標 (子午線に沿う距離)m，xm はより正確
な子午線座標である．
流れの循環 Γ は翼列の外では一定になり，上流側では ΓU = 2πrvuU ，また下流側では Kutta の条件を満足するように予測子修正子法によって決定される．その予測値は

ΓD = 2πrvu = 2πr(wm tan βT +ωr)

(11.50)ただし βT は翼の出口角である．これは翼の出口部における圧力面の方向 βT 0 と負圧面の方向 βT 1 を平均したもので，流速と同様に回転流面上で子午線方向から θ の正の向きに正に測られる．tbt = tan βT は次式によって計算される．
³ β +β ´T0T1tan βT = tan2
周期境界は羽根の前後縁から相対速度の方向に延長される．

(11.51)

η(ξ + ∆ξ)−η(ξ)
wu
=
∆ξ
wm
ただし wu = Γ/(2πr)−ωr である．θ = η ．三角形 1 次要素への分割は，羽根および周期境界の間を周方向に等分割することによって行われる．各節点の上記の番号 ij と左下から振られる通し番号 m の関係は次のようになる．

m = i jp +j +1;

i = (m−1)/jp ,j = MOD(m−1, jp )ただし jp = jf +1，割算は整数演算とする．計算に用いたメッシュをこの節の末尾に示す．
サブルーチン CALPSI： 流れ関数 Ψ の有限要素方程式の係数と右辺を設定し，この方程式を緩和法で解く．まず配列 nm に各要素の節点番号を入力する．なお要素節点番号 1, 2, 3 は左回りに取らなければならな
い．n 番目の要素の nm(1,n)=m とすれば，その要素と上の要素の節点番号は下図のようになる．
m+1 •
@
3

@

@

@

n+1

@
2@
»• m+jl
»
» »3

1 »
»
m• 1

@

@

n

@

@

@2
@• m+jp

30

非粘性流れの FEM 解析

次に下流側の循環 ΓD の修正値を Kutta の条件から求める．下図に示すように，Ψ = 0 の流線はいわゆ
る 滑り のために下方の周期境界よりも通常下側に来る．その原因は偏差角 (deviation) と同じで，流れが羽根によって完全にはガイドされないことによる．この図の T は羽根の後縁点，S, S 0 は 0 流線上の点，

0 は下方の周期境界点とする．Kutta の条件が満足されるときには次の関係が成立する．
1
ηS −ηT = ∆ξ (tan βT +tan βS ) ≡ ∆ξ tan β2
ただし tan βT は既知で，tan βS = wuS /wmS ≈ (wu0 /wm0 )，また付きの量は点 T の量と点 S, S 0 または0 の量の平均を意味する．また下図を参照すれば次の関係が得られる．
(n)

(n−1)

ηS −ηS 0
wu −wu≈ηS −ηT
wu

(n)

(n−1)

ただし wu = (wuT +wu0 )/2，また wu −wu

(n)

(n−1)

は循環の修正値に関係付けられる量で (wu −wu)T ≈ 0とすれば

wu(n) −wu(n−1) =

1 (n)
1
1
(n)
(n−1)
(w −wu(n−1) )0 = (vu(n) −vu(n−1) )0 ≈(Γ −ΓD)
2 u
2
4πrT Dまた ηS 0 は Taylor 展開 Ψ0 = ΨS 0 +(Ψη )S 0 (η0 −ηS 0 ) から

ηS 0 = η0 −Ψ0
wm0のようになる．これらの式から，下流側循環の修正値の式は次のようになる．
(n)

(n−1)

ΓD = ΓD+α∆ΓD ,(n)(n−1)(n−1) ηS −ηS 0∆ΓD ≡ ΓD −ΓD= 4πrT wu
ηS −ηThwm ³
Ψ0 ´i(n−1)= 4πrT wu +ηT −η0 +∆ξ
wm0(11.52)

ただし α はダンピング係数である．

T•

∆ξ
(wm )

AT
AT
ηT −ηS
AT
(wu )
AT
•
0
ηS 0 −ηS
AT
(n)
(n−1)
³
•
³
A
T
)
S 0 ©A•T ³ (wu −wuS
A T Ψ (n−1) = 0
A (n)
AΨ =0

次に有限要素方程式 (11.38) の剛性マトリックスを設定する．剛性マトリックスは式 (11.39) の第 1 式で与えられるが，三角形 1 次要素の場合には式 (11.32) に示すように形状関数の微分は定数，その積分は単に要素面積をかけることになり次のように簡単になる．

φ1,x2 +φ1,y2
φ1,x φ2,x +φ1,y φ2,y∆eAe = [Kkj] = φ2,x φ1,x +φ2,y φ1,yφ2,x2 +φ2,y22
φ3,x φ1,x +φ3,y φ1,y φ3,x φ2,x +φ3,y φ2,y
Ã !
´ φ
∆³
1 T
,x
=
=
Φ Φ
φ ,x φ ,y
2
2∆
φ,y


φ1,x φ3,x +φ1,y φ3,y
φ2,x φ3,x +φ2,y φ3,y 
φ3,x2 +φ3,y2

31

回転流面上の翼列流れのプログラム
ただし ΦT は Φ の転置行列 (transposed matrix) である．あるいは
 eA 11 Ae12 Ae13
Φ11 Φ11 +Φ21 Φ21 Φ11 Φ12 +Φ21 Φ22
1
Ae = Ae21 Ae22 Ae23  =Φ12 Φ11 +Φ22 Φ21 Φ12 Φ12 +Φ22 Φ22
2∆
Φ13 Φ11 +Φ23 Φ21 Φ13 Φ12 +Φ23 Φ22
Ae31 Ae32 Ae33
Φ11 Φ13 +Φ21 Φ23Φ12 Φ13 +Φ22 Φ23 
Φ13 Φ13 +Φ23 Φ23(11.53)

剛性マトリックス A は，図 11.2 に示すように帯行列になり太線の外の要素はすべて 0 になる．ここではA に記憶することにする．その配列 A への入力は要素ごとに行われる．
A のバンド幅の中だけを縮小行列 Ã
例えば n = 81 番目の要素の要素剛性マトリックスは，図に示す三角形要素の 3 節点の値を用いて計算され，A に既にある値に重ねて下表に示すように入力される．すなわちこの入力は図に ◦ 印得られた 9 つの値は Ã
を付けて示す 9 要素に対して行われることになる．一方 連立 1 次方程式の例えば 40 番目の式は，形状関A の 7 要素にのみ値が入り，残りは 0 のままである．
数 φ40 を重み関数とする式で，その係数は下表に示す Ã
その値は n = 67, 68, 69, 80, 81, 82 の 6 要素の計算で得られた要素剛性マトリックスの値の和になる．
1

j

1

−lf 0
1

K

@ @
@ @
@ @@ @◦
40 @ ◦ ◦@@ ◦@@◦ ◦ @
47
◦ ◦◦
48
@
@ @
@ @ @
@ @@ @
@ @ @
k
@ @ @
@ @
@ @
K

l

lf

41
»•
33 •»» @@ 69 @@ 82 »»@• 4868 @40•»»»»32 •@ 67 @ 81
@
80 @
@
•
»
@•» » 47
48
39
»»•

¡
¡

◦ ◦◦
40
¡
◦¡
◦
◦¡
47 ¡
◦¡
48 ◦ ◦¡

k=k
l = j −k
⇐⇒

40 •»

@ 81
@
@•

k

¡
¡

K

A

A
Ã

47

A
図 11.2: 三角形 1 次要素の剛性マトリックス A とその縮小行列 Ã
l = −8−7−1

m = 4047
48

n = 67
68
69
80
81
82

Ae31
Ã40,−8
Ae31
Ae21

Ae21

Ã40,−7
Ae23
Ae21

Ae32

0
Ae11

1

Ae22
Ae33

Ae23

Ã40,−1
Ae32
Ae31

Ã40,0
Ae33
Ae22
Ae22
Ae33
Ae11
Ae11

7
Ae12

Ã40,1
Ae23
Ae13

8
Ae13

Ã40,7

Ae32
Ae12

Ã40,8

Ae13
Ae12

有限要素方程式 (11.38) の右辺の外力項は，その係数，式 (11.39) の第 2 式を三角形 1 次要素の積分公式(11.33) を用いて計算することによって設定される．すなわち式 (11.40) の右辺 be は次のようになる．

2f1 +f2 +f3∆e
be = −(Fkjfj ) = − f1 +2f2 +f3 (11.54)24
f1 +f2 +2f3
ただし f は流れ関数の方程式 (11.21) の右辺で，この羽根車では式 (11.23) で与えられる．上記のプログラ
ムではこの右辺 B の計算を配列 A の計算と同時に行っている．なお質量の集中化 (mass lumping) は行わ
ない．

32

非粘性流れの FEM 解析

ここで境界条件の処理についてまとめて説明する．羽根表面上の節点の有限要素方程式は Dirichlet 条件の式 (11.25a) とし，Ãml = 1 (l = 0), = 0 (l 6= 0), bm = 0 or q と置く．上流と下流境界には Neumann 条
件 (11.25b) が課される．これらの境界上の節点の有限要素方程式 (11.38) は，内点とは異なるものになり，その剛性マトリックス Ãml と外力項 bm の値は境界の片側にある要素からのみ入力され，また右辺には自然境界条件が加わる．この境界条件の係数は式 (11.39) の第 3 式で与えられるが，三角形 1 次要素では次のようになる．

(∆s/3 (j = k)φk φj ds =e∆s/6 (j =
6 k)
∂Ω2

Z

(11.55)

ただし ∆s は境界要素の長さである．また Ψ の法線微分値は上流境界では wu ，下流境界では −wu となる．
プログラムではこれらの境界点の有限要素方程式の右辺に ±wu ∆s を加えている．周期境界では，まず上方境
界点の有限要素方程式を下方境界点の Ãml , bm と周期境界条件 (11.25c) を用い完全なものに作り変える．また下方境界点の有限要素方程式は周期境界条件の式とし，Ãml = 1 (l = 0), = −1 (l = jf ), = 0 (otherwise),bm = −q と置く．
Ψm の連立 1 次方程式は 14.1.4 項に述べた緩和法で解かれる．緩和計算は i1 ≤ i ≤ i2 の翼列領域，
0 ≤ i < i1 の上流領域，i2 < i ≤ if の下流領域に別けて行われる．翼列領域では，上流から下流へ，下流から上流へ向かう線過緩和の掃引が交互に行われる．上流領域では，まず式 (11.44) によって列ごとに一様な修正値が加えられ，それから上流へ向かって線過緩和の掃引が行われる．ただし上流境界は線緩和しない．
下流領域の緩和計算も同様に行われる．最後に緩和計算の収束状況を調べるために，有限要素方程式の残差を計算し，最大最小残差値を求める．
サブルーチン VELCTY： 流れ関数 Ψ を微分し流速成分 wm , wu を求める．その方法は 14.1.5 項に述べたもので，各節点まわりに Ψ (ξ, η) の近似 2 次多項式を重み付き最小２乗法によって作り，この式を微分し流速成分を求めるというものである．その際に必要な 6 近傍節点の節点番号は，あらかじめサブルーチン
SETNEB で配列 NEB に入れておく．ここでは各節点ごとに，式 (11.46) の pb(k) = Ψk − Ψ0 と ps(l,k) = ψlk ，式 (11.48) の a(n,l) = anl と f(n) = fn ，逆行列 a(l,k) = (a−1 )ln の計算が順に行われる．ターボ機械羽根車の非粘性相対流れでは，岐点エンタルピーに相当するロータルピー rot = I = p/ρ+v 2/2−ωrvu は流れ
場を通して一定になる．したがって静圧は次式から求められる．
pv2(ωr)2 w2
= I − +ωrvu = I +
−
ρ
2
2
2
ただし I = p0U/ρ−ωΓU/(2π) である．
６つの補助的サブルーチン： ESPL は，指数スプラインで補間するサブルーチンで，配列 xs,us にスプ
ライン曲線 u(x) 上の基準点の x と u の値，ls に基準点の番号 i，g1,gf にスプライン曲線の両端の勾配

u0 (x)，x に基準点を含む補間点の xi を入力しこのサブルーチンを引用すれば，配列 u に補間値 ui が得られる．GAUSS3 は３重対角行列の連立 1 次方程式を解くサブルーチンで，a に縮小行列，b に右辺を入れて引用
すれば，b に解が得られる．また GAUSSP は GAUSS3 の周期条件が課された場合の version である．SETNEB
は，各節点まわりの 6 近傍点の節点番号を配列 neb に入力するものである．WRTF は回転流面上の節点値を出力するサブルーチンである．例えば静圧の計算値を 10 倍したものを出力したい場合には，メインプログ
ラムで WRTF を，関数 f = p，表題 z1 = p*10 ，桁移動数 isf = 1 とおいて CALL する．なお桁移動数は−3 から 3 まで取ることができる．MATINV は配列 a の行列の逆行列を同じ配列 a に出力するサブルーチン
で，ピボットの絶対値が e よりも小のときには行の入替えが行われるが，e よりも大にできないときには計算がストップし 5555 が出力されるようになっている．

回転流面上の翼列流れのプログラム

33

反復計算は下表に示すように，n = 60 でほぼ完全に収束するが，実用上は 20 回程度でも十分であろう．
なお有限要素方程式の最大最小残差値 Resmax , Resmin は，解が収束に近付いた時点では，下流境界点を除いてほぼ 0 になる．図 11.3 に計算に用いたメッシュと，計算で求めた圧力分布と流線を示す．
n
1
2
10
20
40
60

Γd
2.803
2.651
2.100
1.900
1.799
1.791

Resmax
.0197
140
60
22
21
20

Resmin−.0277
−178
−51
−24
−26
−25

図 11.3: ２次元翼列流れの計算メッシュ，圧力分布と流線

34

非粘性流れの FEM 解析

11.4

３次元翼列ポテンシャル流れの FEM 解析

飛行機の翼まわりの流れには後流渦による吹き降ろし (down-wash) があり，そのために翼の揚力が低下
する．飛行機周りのポテンシャル流れ解析では，通常翼の下流側に随伴渦面 (trailing vortex sheet) が考慮される．ターボ機械の羽根車を通る３次元ポテンシャル流れ解析でも，翼負荷分布一定の自由渦設計の場合を除いて，随伴渦面が考慮されるべきである．以下に述べる３次元翼列ポテンシャル流れの解析は，水車ランナー11 や亜音速羽根車12 の FEM 解析に随伴渦面を考慮し，下流側の流れを合理的に計算しようとするものである．

11.4.1

ポテンシャル方程式の境界値問題

羽根車の絶対流れは，上流の静止流体に由来する渦無し (irrotational) 流れで，また相対流れは定常と仮定できる．また完全気体の等エントロピー流れを仮定すれば，この流れの基礎方程式は次のようになる．

∇·ρ w = 011
w = ∇w2 −w
w ×∇×w
w = − ∇p+rr ω 2 −2 ω ×w
w
w ·∇w
2
ρ
T ds = 0 = dh−dp/ρ

(11.56a)(11.56b)(11.56c)

ただし v = w + ω ×rr で v , w , r , ω はそれぞれ絶対速度，相対速度，半径，羽根車の回転の角速度である．
また ρ, p, T, s, h は流体の密度，静圧，絶対温度，比エントロピー，比エンタルピーである．非圧縮性流れの場合には，これらの式において ρ = 1，また第 3 式は不要になる．
式 (11.56b) を相対流れの 1 つの流線に沿って積分し13 ，速度線図から得られる関係

w2 = v 2 +(ωr)2 −2ωrvu
を用い，また圧縮性流れの場合には式 (11.56c) を用いれば，次の関係が得られる．

p v2


 + −ωrvu(非圧縮性流れ)I ≡ ρ 22
 h+ v −ωrv

(圧縮性流れ)u2(11.57)(11.58)

これより上流の全圧 P = p+ρv 2/2 または岐点エンタルピー H = h+v 2/2 が一様ならば，上式で定義される
ロータルピー (rothalpy) I は流れ場を通して一定になる．
絶対流れは渦なしであるから速度ポテンシャル φ を次のように定義することができる．

v = −∇φ(11.59)また圧縮性流れでは，式 (11.56c)，音速 c の式

c2 = (∂p/∂ρ)s = (k − 1)h(11.60)

11 T. Nagafuji and H. Morii, Proc. IAHR 10th Symp., Vol.1(1980), 583.
12 T.E. Laskaris, AIAA J., 16(1978), 717.
13 流線方向の単位ベクトル 1 とのスカラー積を作り流線に沿って積分する．同じ方向の 2 つのベクトルを含むスカラー３重積は 0
RRRw ) = 0，1 · ω ×w
w = 0．また 01 ·∇a ds = a−a0 ， 01 · r ds = 0 r dr = (r2 −r02 )/2．
であるから，1 · w ×(∇×w

３次元翼列ポテンシャル流れの FEM 解析

35

式 (11.56a) を順に用いれば，次の関係が得られる14 ．

11
w = −c2 ∇·vv
w ·∇h = w ·∇p = c2w ·∇ρ = −c2 ∇·wρρ
式 (11.59) にこの関係を用いれば，次の速度ポテンシャルの方程式が導かれる．
(
0(非圧縮性流れ)2∇ φ=b≡
2w ·∇h/c(亜音速流れ)(11.61)

この方程式は基礎方程式 (11.56) の代わりをなすもので境界値問題として解かれる．流れ場の φ が決定できれば，これを数値微分することによって速度 v ，w ，ロータルピー一定の条件 (11.58) からエンタルピー h，
式 (11.60) から音速 c, . . . を順に求めることができる．つまりこの境界値問題は，亜音速流れの場合には，
式 (11.61) の右辺に適当な予測値を与え予測子修正子法によって解かれることになる．
計算領域は１つの翼間流路とその上流境界と下流境界までの延長部分からなる．羽根車と流路の形状，羽根車の回転の角速度のほかに，上流と下流の条件を過不足なく与えることが肝要である．例えばポンプまたは送風機を考えれば，流量と予旋回 (prerotation，流入流れの旋回) を与えることができる．ここでは上流
境界の子午線速度 wmU と循環 ΓU を与える．なお絶対流れは渦なしであるから，上流領域では Γ = 2πrvu
一定の自然渦になる．更に圧縮性流れでは熱力学的状態量 (properties)ρ, p, T, s, h, . . . またはそれに代わる c, M, h0 , . . . の中から 2 つを選び与えなければならない．ここでは上流境界の岐点密度 ρ0U と相対 Mach
数 MrU を与える．wmU と ΓU が与えられれば残りの速度成分は速度線図から求めることができる．またρ0U と MrU が与えられれば，残りの熱力学的状態量は次の完全気体の状態方程式，Mach 数の定義式，等
エントロピー流れの式などから決定することができる．
p = ρRT = (k−1)ρh/k,

h = Cp T,² = Cv T

M = v/c, Mr = w/c³ h ´1/(k−1) ³
ρ0
k−1 2´1/(k−1)
0
=
= 1+
M
ρ
h
2
ただし R, Cp , Cv , k = Cp /Cv は気体定数，定圧比熱，定積比熱，比熱比，また ² は比内部エネルギーである．下流側のこれらに相当する量はすべて計算によって求められるが，上流の条件の一部を下流の条件で換えることももちろん可能である．
式 (11.61) の境界条件は次のように与えられる．

Dirichlet 条件： ∂Ω1
φ(z, r, θ) = φ(z0 , r0 , θ0 ) − ΓU (θ−θ0 )/2π

(上流境界)(11.62a)

Neumann 条件： ∂Ω2



 −wm
φn ≡ n ·∇φ = f ≡ ±nn · ω ×rr


0

(下流境界)(羽根表面)(11.62b)(ハブ・シュラウド・ケーシング面)

周期境界条件： ∂Ω3

φ(z, θJ , r) = φ(z, θ0 , r) − Γ/Nφn (z, θJ , r) = −φn (z, θ0 , r)

(上流・下流領域の周期境界)(11.62c)

14 多少蛇足を加えれば，熱力学の関係式は流体の系 (system) に対して成立するもので，その微分演算子は絶対流れに対しては
˜
d/dt ≡ ∂/∂t+vv ·∇，また定常の相対流れに対しては d/dt
≡ w ·∇ のように置くべきものである．
˜v = w ·∇ρ+ρ∇· w = 0，しかるに dρ/dt = dρ/dt
また連続の式は dρ/dt+ρ∇·vであるから ∇· v = ∇· w である．これらの式は同じ現象，流跡線に沿う密度変化と流れの発散を異なる座標系から見ているものである．なお ∇· ω ×rr = 0 である．

36

非粘性流れの FEM 解析

ただし ∂Ω = ∂Ω1 +∂Ω2 +∂Ω3 は計算領域 Ω の境界，θJ = θ0 +2π/N ，N は羽根枚数，下添字 n は境界における法線微分 ∂/∂n = n ·∇，n は外向き単位法線ベクトルである．
上流または下流の 流れが周方向に一様になるところでは次の子午線流速 wm と循環 Γ の関係が成立する．
∂wm
1 n ω ∂Γ
1 ∂Γ 2 o=−(11.63)∂s
2wm π ∂s (2πr)2 ∂s
また流量と質量流量は次式で与えられる．
Z
Z
Q = 2π wm r ds = const.
G = 2π ρwm r ds = const.

(11.64)

ただし m, u, s は子午線方向，周方向，それらと直交する方向の長さまたは成分を意味する．次に式 (11.63)の導出過程を示す．式 (11.56b) の s 方向成分の式で vs = 0 と置けば，一般化された径方向均衡条件 (radialequilibrium condition) の式が得られる．
1
1 ∂p= vu2 cos γρ ∂sr
ただし γ = tan−1 (dr/dz) は子午面流線の回転軸に対する傾きである．この式にエネルギー方程式 (11.56c)とロータルピー一定の条件 (11.58) を用いれば

vm

∂vm
∂rvu
1
∂vu
=ω
− vu2 cos γ − vu
∂s
∂s
r
∂s

さらに Γ の定義式を用いれば上式が導かれる15 ．なおこの式は actuator disc 理論の等エントロピー流れ
の式

o
dP
1n
1 d
= 2 −rζu +(rvu )2dψr
2 dψ(11.65)に相当のものである．上流境界では絶対流れは渦なしで，ΓU 一定であるから，式 (11.63) より wmU も一定
で，流入速度は一様になる．他方下流境界では一般に後流渦面が存在するので ΓD は一定にならず，ΓD の
値は Kutta の条件から決定され，wmD の値は式 (11.63) と (11.64) を共に満足するように計算によって求められる．ここでは式 (11.61) の境界値問題は，六面体３重 1 次要素を用い変分 FEM によって解かれる．

11.4.2六面体３重 1 次要素の補間関数

ここでは x 空間内の六面体要素を ξ 空間内の立方体要素 (−1 ≤ ξ, η, ζ ≤ 1) に写像し，すべての計算を写
像空間の立方体要素に関して行うことにする．立方体要素の補間関数は ξ, η または ζ に関して 1 次で，その形状関数は値が当該節点で 1，他の節点で 0 でなければならないから次のようになる．

ψk (ξ, η, ζ) =1
(1+ξk ξ)(1+ηk η)(1+ζk ζ)
8(11.66)

ただし ξk , ηk , ζk の値は，要素節点番号を下図に示すように振れば下表のようになる．
k

1

2

3

4

5

6

7

8

ξk
ηk
ζk

−1
−1
−1

1
−1
−1

−1
1
−1

−1
−1
1

−1
1
1

1
−1
1

1
1
−1

1
1
1

15 v cos γ +r∂v /∂s = ∂(rv )/∂s．
uuu

３次元翼列ポテンシャル流れの FEM 解析

37

5

»»»HH
4 »»»
HH
H
HH ζ
»8
»H
HH»6
»»»
6

:η


HH
j
»»»3HHξ
»
»
»
HH
k=1 H
HH
»7
»H
»
HH»»»
2

この形状関数の ξ, η, ζ に関する微分は次のようになる．

1
ξk (1+ηk η)(1+ζk ζ),
81
ψk,η = (1+ξk ξ)ηk (1+ζk ζ),
81
ψk,ζ = (1+ξk ξ)(1+ηk η)ζk
8
ψk,ξ =(k = 1, 2, . . . , 8)

これを具体的に書けば
 

ψ ,ξ
−(1−η)(1−ζ) (1−η)(1−ζ)
  1
ψ ,η  = −(1−ζ)(1−ξ) −(1−ζ)(1+ξ)
8ψ ,ζ
−(1−ξ)(1−η) −(1+ξ)(1−η)··················(11.67)(1+η)(1+ζ)
(1+ζ)(1+ξ)
(1+ξ)(1+η)

また形状関数の x, y, z に関する微分は次式から求められる．

 
 
ψ ,ξ
ψ ,x
 
 
ψ ,η  = Jψ ,y 
ψ ,ζ
ψ ,z

(11.68)

ただし J はヤコビ行列でアイソパラメトリック要素では次のようになる．

x1 y1

x2 y2

x y

  
 3
3
xξ yξ zξψ ,ξ ³
´
³
´ 

x
y

  
4
4
J = xη yη zη  = ψ ,η  x y z ,x y z =x y
 55xζ yζ zζψ ,ζ
x6 y6

x y
7
 7
x8 y8


z1

z2 

z3 


z4 

z5 


z6 

z7 

z8

(11.69)形状関数の 1 つの六面体要素にわたる積分は，立方体要素にわたる積分として計算される．例えばZ

Z 1Z 1Z 1x)ψl (xx)dΩ =ψm (xΩeψm (ξξ )ψl (ξξ )J(ξξ )dξdηdζ =−1 −1 −1

ただし J = |J| はヤコビアン，ξ k は Gauss 点の座標である．

8
X

J(ξξ k )ψm (ξξ k )ψl (ξξ k )k=1(11.70)38

非粘性流れの FEM 解析

11.4.3

変分 FEM によるポアソン方程式の解法

ポアソン方程式 (11.61) と境界条件式 (11.62) の境界値問題はここでは変分 FEM で解かれる．この境界値問題に相当の変分問題の汎関数は，式 (11.61) の b と式 (11.62b) の f の値が当面 確定しているものとす
れば次のようになる．
Z n
Z
o
1
J(φ) =(∇φ)2 +bφ dΩ −f φ dSΩ 2∂Ω2(11.71)

変分 FEM では計算領域を要素に分割し，各要素内で φ, b, f を次の補間関数で近似する．

x) = ψm (xx) φm ,
φ(xx) = ψm (xx) bm ,b(xx) = ψm (xx ) fmf (x(11.72)x) は形状関数，φm , bm , fm は節点値である．また積分は要素ごとに実行される．このとき汎関
ただし ψm (x数の式 (11.71) は次のようになる．

J(Φ) =

N ³
X
1
n=1

2

e
e
Hml
φm φl + Kml
φm bl − Geml φm fl´(11.73)

ただし

Z
e
Hml
=

Z
ψm,i ψl,i dΩ,

Ωe

e
Kml
=

Z
ψm ψl dΩ,ΩeGeml =ψm ψl dS

(11.74)∂Ω2eRRu,i = ∂u/∂xi ， Ω e dΩ は要素ごとの積分， Ω e dS は Neumann 条件の課されている境界の要素ごとの積分，ま2
P
た はすべての要素の和を意味する．汎関数 J をすべての節点値に関して最小化するために ∂J(Φ)/∂φm = 0と置けば次の有限要素方程式が導かれる．
N
X
n=1

e
Hml
φl = −

N
X
n=1

e
Kml
bl +

N
X

Geml fl(m = 1, 2, . . . , M )(11.75)n=1

ここでは六面体３重 1 次要素 (hexahedral trilinear element) を用いることにする．式 (11.74) の係数の値
は，六面体要素を立方体要素に局所写像し，この立方体要素上で見積もられる．式 (11.75) に上流境界 ∂Ω1の条件 (11.62a) と周期境界 ∂Ω3 の条件 (11.62c) を考慮したものは正則な連立 1 次方程式になり，節点値

φm の値はこの方程式を解くことによって一意的に決定される．ここでは反復法によらず，ガウス消去法で直接解くことにする．なおその際に，桁落ちを防ぐために必要に応じ，倍精度実数型でプログラミングする
こと，流量が下流境界まで一定に保たれるように wm を補正すること，更にゆがみの小さい要素を用いることも検討しなければならない．次に随伴渦面の考慮の仕方，精度の良い数値微分，上流と下流域に２次元解析を補う３つの改良点について述べる．

11.4.4随伴渦面を考慮した下流領域の計算

随伴渦面 (trailing vortex sheet，後流渦面) は羽根の後縁から脱落する渦糸 (vortex filaments) から構成
され，相対流れに沿って形成される．この渦面の両側で相対流速は不連続になるが，ロータルピー一様の流
れでは方向のみ異なり大きさは等しくなる．また相対流れは定常で渦糸は渦面両側の相対流速の平均速度の方向に伸びることになる．計算を簡単にするために，周期境界をこの渦面上に取り，下流領域を周期境界と渦糸を含む回転面 (ケーシング) に合わせて要素に分割することにする．下流領域の循環 ΓD は各回転面

３次元翼列ポテンシャル流れの FEM 解析

39

上で一定になり，その値は Kutta の条件から決定される．下流領域の要素分割と，循環 ΓD (s) と子午線流
速 wmD (s) の計算は予測子修正子法によって行われる．
予測子の計算： 回転面 (k 面) は流路幅を等分割する面に取られ，ΓD は各回転面上で一定とし式 (11.76a)
から求められ，wmD は下流境界で一様とし式 (11.76b) から求められる．周期境界 (j0 面) は翼後縁から式 (11.76c) によって延長される．

ΓD = 2πr(ωr+wm tan βT )Z
G = 2πρU wm r ds

(11.76a)(11.76b)´
r dθ
wu
1 ³ ΓD
=
=
−ωr
dm
wm
wm 2πr

(11.76c)R
ただし βT は翼後縁における回転面上の反り曲線の方向で子午線方向から測られたもの，2π r ds は流路断面積である．
修正子の計算： 各回転面は式 (11.77a) により１つの渦糸を含むように取られ，ΓD は式 (11.77b) により

Kutta の条件を満足するように修正される．なおこの修正子にはダンピングをかけることが必要である．
wmD は式 (11.63)，周期境界は式 (11.77c) によって修正される．
wr
dr
=
dz
wz

(11.77a)δΓD = 2πr(wm tan βT −wu ),

(n+1)

ΓD

(n)

= ΓD + αδΓD

r dθ
wu
=
dm
wm
ただし 上添字(11.77b)(11.77c)は渦面上の平均相対流速を意味する．収束を加速するためには ΓD だけでなく，回転面と周期境界の計算でも適切なダンピングをかけることが望ましい．収束解は通常 5 回程度の反復で得られる
が，高マッハ数の亜音速流れでは更に多くの反復が必要である．
全体の計算手順は次のフローチャートのようになる．

11.4.5

FEM における精度の良い数値微分―３次元

Laplace 方程式の 8 節点立方体要素を用いた FEM による３次元解析は，2 次中心差分のものよりも高精度であることが知られている．しかしながら一般の六面体要素では，直交性や滑らかさの問題のために精度が低下し，極端に歪んだメッシュでは計算が不能になることすらある．一方，関数の微分値は，3 点差分
式では 2 次精度であるのに対し，六面体要素の FEM では各要素内で補間関数を微分したのでは 1 次精度に留まる．ここには，8 節点六面体要素または 6 節点三角柱要素の FEM 解析で微分値をより良い精度で求める方法について述べる．
この方法は，すでに 11.1.4 項に述べたものと基本的に同じで，その式 (11.45)–(11.48) は３次元の場合にもそのまま成立する．ただし 式 (11.45) で，L = 9, ψ1 = x1 , ψ2 = x2 , ψ3 = x3 , ψ4 = x12 , ψ5 = x22 , ψ6 =x0 は節点 x 0 に原点を持つ局所座標である．また 9 つ
x32 , ψ7 = 2 x2 x3 , ψ8 = 2 x3 x1 , ψ9 = 2 x1 x2 , x = x −xの未定係数 cl は，点 x 0 近傍の K = 14 個の節点条件 (11.46) が平均的に満足されるように，最小２乗法によって決定される．その方法は前項に述べたものと同じになるのでここでは省略する．重み関数も前項同様に Wk 2 = ( x 2 )−3 と置く．14 近傍点は図 11.5 に示すように，すべての方向に均等に 3 節点ずつ並ぶように取られる．なお近傍点の 14 の番号は，特に隅の点で矛盾を来たさないよう十分考えて振らなければならな
い．この方法によれば微分値を，良質のメッシュに対し 2 次相当精度で求めることができる．

40

非粘性流れの FEM 解析?

¨¥§START ¦Input
Cascade Geometry
and Flow Conds

Compute b- n = n+1

n=0

Correct k-surfaces,ΓD ,wmD
and Vortex Sheet

Discretize Mesh

Compute φ

Determine Upstream
Bound Conds andPredict Downstream
Bound Conds

Compute w , M, p¢¢ WRITE ¢¢¡Does Solution @
¡
@
@
¡ No
@Converge?
¡
Yes¢¢ WRITE ¢¢¨¥§STOP ¦
図 11.4: 3 次元翼列流れ FEM 解析のフローチャート

11.4.6上流と下流域に２次元解析を補う方法

普通のポンプの羽根車は翼列の食違い角が非常に大きいため，子午面上の格子線の回転面 (i 面と k 面) と翼間を分割する面 (j 面) から普通に作られた六面体要素は，ゆがみの非常に大きいものになる．三角柱要素を用いれば，ゆがみをある程度小さくできるが，それでも羽根車から上流または下流に離れるとゆがみは大きくなり，解の精度低下，更には計算不能に陥ることもある．この項には羽根車から離れた上流と下流領域に軸対称２次元解析を補う方法について述べる．この方法では計算領域を羽根車を含む３次元流れ解
析領域と，その上流と下流の２次元の軸対称流れ解析領域に分け，３領域を交互に計算する．この方法ではプログラムは多少煩雑になるが，上流側と下流側の十分広い領域を計算に含めることができ，また中心軸上の流れも問題なく計算できるという利点もある．
この２次元軸対称流れの支配方程式は次のように書くことができる．

∂ ³ 1 ∂ψ ´ ∂ ³ 1 ∂ψ ´+
= ζu∂z rρ ∂z∂r rρ ∂r(11.78)ただし ψ はストークスの流れ関数で

vz =

1 ∂ψ
,
rρ ∂r

vr = −1 ∂ψrρ ∂z(11.79)のように定義される．また ζu は絶対流れの渦度の周方向成分で，上流領域ではゼロ，下流領域では随伴渦に相当する周方向に一様な渦度で次のようになる．

ζu = −ωrd ΓD 1 d ³ ΓD ´2+dψ 2π 2r dψ 2π(11.80)

なおこの式はロータルピー一定の条件 (11.58) のもとでは actuator disc の理論の式 (11.65) と等価である．

41

３次元翼列ポテンシャル流れのプログラム
14
»»•H
»
H
»
H •6»»»HH•
»» H
12 •H
» 10
H »»»HH
HH 8 »
»»»
»
»
»
•
H
4•»
H •2
»»H
H
»»»H
»
◦0» 3 »H
H
H
»
»
HH
»
H
•»
»
k = 1• H»»
»
•H
»
» 13H
H •5»»»HH•
»»»H
11 •H
»
HH
» »HH
»»» 9
H•»»»

10
•H

14

11
•H

13
5»»»•
H
•
»
9 »
10 •»
14
•H
»»•HH
H
»
»
H
12
»
»»
H
•»
•6»
»
HH»»
»
•H
»
» 8H
7 »»•

H •4»»»HH
3H
•»

•
H
»
HH
»»»H
»
»
H
2
»
◦
» 0
H•»

k=1

H »»»•HH
»»
»»»HH
12 •H
•6»
»
•H
»»
H»
»»
3
8H
•
H •4»»»HH•2
»»»H
7 •H
H »»»H
H 0»»»
HH»
1 »»◦
•
13
9 •H
H
»»»•HH»H»»»
H•5» »
11 •H
H»»»

7
(a) 六面体要素の通常点(b) 境界面上の点
• 13

¡
(¡
P
((

(PP
12(((
P• 14
•P 11¡
(((
P
(
•
(
¡
P
¡
(
¡
(
¡
P
P
¡
(
PP
(
13(
P
•¡(
6
k = 1H
(
•
( PP
H•(
P• 5
P(
¡ ((
P0(
◦
(
¡
P(
¡
(
¡
¡
P
P
¡
2 •(
(
P
•
(
P
4
P•(
¡(
3 (
10
(
•
P
(
( ¡ (P
P
(
P
P(
•P(
(
¡
¡
(
¡
¡
(
P
P
¡
8 •PP
(•
7 (
P(
¡(
9
(a) 三角柱要素の通常点

P
P• 14
12 •(
(
PP (
¡((
•P
¡
P¡
11 P
•3
¡
6
(
•
(¡
PP
((• 2
P(
1 •(
P(
•5(
P◦P
(((
¡
¡
P•¡
0 P
4
10
(
•
(
P
(
(
P
P((((• 8
PP (
¡((
¡
•P
¡
P•¡
7 P
9
(b) 境界面上の点(c) 隅の点
10(•PP (((• 8
(((
P(
P
¡((
P•P
(
¡
¡
P9•¡
7 P
©13
•©
¡
(¡
PP
(((
P
12 •P
(• 14
P•P
(((
¡
¡
P¡
11 P
•3
¡
6¡
P
(
PP (((•P 5(((• 2
P
P(
•P
PP ¡
•
((
(
k=1 P
◦P
P
¡
P
P
¡
P
•P
0
P
4 P
(c) 隅の点

図 11.5: 14 近傍節点の選び方
式 (11.78) の境界条件は次のように与えられる．

Dirichlet 条件： ψ = 0

または

G/2π

Neumann 条件： ∂ψ/∂n = ∓rρvs
接合条件：
ただし 上添字±1 ∂ψ
∂φ= vs = ws = −
rρ ∂n
∂s

(ハブ・ケーシング面)(11.81a)(上流と下流の境界)(11.81b)(３次元流れ領域との境界)(11.81c)は３次元流れの周方向平均値である．流れ関数 ψ の方程式 (11.78) と境界条件 (11.81) の境界値問題は，四辺形双 1 次要素を用い，変分 FEM で容易に解くことができる．なお３次元流れ領域の計算における軸対称流れ領域との接合条件は次のようになる．

∂φ1 ∂ψ= −v n = −wn = ±
∂n
rρ ∂s

(11.82)全体の計算は，上流の軸対称流れ → ３次元流れ → 下流の軸対称流れの順に反復される16 ．

11.5

３次元翼列ポテンシャル流れのプログラム

メインプログラムは，データの入力，サブルーチンの引用，計算結果の出力部分からなる．このプログラ
ムでは，円柱座標系 zθr とデカルト座標系 x1 x2 x3 が用いられる．x1 = z, x2 = r sin θ, x3 = r cos θ．子午面上の流路の形状 r = r(z) と翼面の形状 θ = θ(z, r) は zθr 座標を用い離散的に与えられ，FEM の計算はす
べて x1 x2 x3 座標で行われる．FEM の計算には 8 節点六面体要素の格子状のメッシュが用いられ，節点には2 系統の番号 ijk と m が振られる．i は上流から下流，j は翼間，k は翼スパン方向の番号，m は通し番号
16 大宮司・高野，上流と下流の効果を含めた三次元ポテンシャル翼列流れの有限要素法解析，日本機械学会論文集，52(1986),
3139–45.

42

非粘性流れの FEM 解析である．m = (i kp +k) jp +j+1，逆に i = (m−1)/(jp kp ), j = MOD(m−1, jp ), k = MOD((m−1)/jp , kp )．
ただし jp = jf + 1, kp = kf + 1，割算は整数演算とする．i = 0：上流境界，i1 = 6：翼前縁，i2 = 14：翼後縁，if = 20：下流境界，j = 0：翼の圧力面，jf = 6：負圧面，k = 0：ハブ面，kf = 6：ケーシング面．
羽根車の回転角速度 ω = 1.0，羽枚数 nblade = 6，上流境界の循環 ΓU = 0，子午線速度 wmU = 0.566，
相対マッハ数 MrU = 0.7，岐点密度 ρ0U = 1.0 を与える．比熱比 k = 1.4．諸量は適宜無次元化されているものとする．計算の手順は概ね図 11.4 のフローチャートと同じである．計算の反復回数 n = 10，解の収束は子午面上の渦糸の r，周期境界の θ，循環 ΓD および流れ場の相対マッハ数 Mr の最大変化量
(n)

(n−1)

Drmax,Dthmax,DGmax,DMmax から判定される．Drmax = max |rik −rik|, . . . ．最後に，計算で求めた速度ポテンシャル φ，相対速度 wz , wu , wr ，静圧 p，相対マッハ数 Mr が出力される．なおプログラムは本節の終りの方にある．
サブルーチン SETX： ここでは計算領域を六面体要素に分割する．まず子午面上の流路が，ハブとケーシ
ング上の与えられた点を結ぶスパン方向の線と，この線を等分割する点を結ぶ流線方向の線によって四辺形要素に分割される．子午面上の流路幅 br = b は，子午線方向の距離を m，スパン方向の距離を s とすれば，近似的に b = nz sz +nr sr で与えられる．ただし (sz sr ) は子午面上の翼スパンの長さベクトル，(nz nr ) = (−∆r ∆z)/∆m は子午線の単位法線ベクトルである．なお上流と下流の境界は流路中央子午線に対し直交するように取られている．
次に翼面が 2 次補間によって形成される．それからこの翼面に繋がる周期境界が次の手順で形成される．
翼の出口角 βT は翼後縁部における圧力面の方向 βT 0 と負圧面の方向 βT 1 の平均とし，子午面上で子午線方向から測られ θ の正方向に正とする．tbt = tan βT は次式によって計算される．
tan βT = tan³β

T 0 +βT 1´

2

ただし βT 0 と βT 1 は翼後縁の丸み rt = rT = .005 を考慮し次のように見積もられる．

βT 0 = tan−1(r∆θ)0 +rT cos βT,∆m0 −rT sin βT

βT 1 = tan−1(r∆θ)1 −rT cos βT∆m1 +rT sin βT

下流側の流れの循環 Gmd = ΓD は次式から求められる．

ΓD = 2πrvu = 2πr(wu +ωr),
wu = wm tan βT ,wm = c(s)wmU

rU bU
rb

(11.83)

ただし c(s) は係数で，翼負荷分布が自由渦設計の場合には 1 であるが，半強制渦設計の場合には子午線速度がハブ側で小，ケーシング側で大になるように取られる．
上流側の周期境界は翼前縁から相対速度の方向に次式により延長される．

wur{θ(m + ∆m)−θ(m)}=∆mwm
下流側の周期境界も同様に翼後縁から相対速度の方向に延長される．六面体要素への分割は翼間および周期境界間の流路を周方向に等分割することによって完了する．なお COMMON 文のブロック名は，メインプロ
グラムに関わるものを除き，その配列や変数の値が決定されるサブルーチン名に因んで付けられている．
サブルーチン PREDCT： 上流境界では諸量は周方向に一様とする．その流路中央における値は与えられた

３次元翼列ポテンシャル流れのプログラム43

データ ΓU , wmU , MrU , ρ0U をもとに次式によって求められる．

p
ΓU
w
,
wu = vu −ωr,w = wmU 2 +wu2 ,
c=
,
2πr
MrU
c2
wmU 2 vu2
w2
(ωr)2
h=
,
h0 = h +
+
,
I =h+
−
,
k−1
2
2
2
2
³ h ´1/(k−1)
X
k−1
p0 =
ρ0U h0 ,ρ = ρ0U,
GU = 2πρwmU r∆skh0vu =なお 3 行めの第 2 式は等エントロピー流れの式である．また上流境界の速度ポテンシャルと相対速度は次式によって与えられる．

Γ(θ−θ∗ ),2π
Γubrwz = wmU ,
wu =
− ωr,
b
2πr

φ(θ) = 1 −wr = −wmUbzb

ただし bz , br は流路幅 b の z, r 方向の長さ成分である．下流境界の子午線流速 wmD は，式 (11.83) から求められた後，下流境界の流量が上流境界の流量 GU に一致するように補正される．
サブルーチン CORRCT： 下流領域の修正計算を回転面 k 面，循環 ΓD ，周期境界の θ，下流境界子午線速度wmD の順に行う．k 面の修正は，この面が 1 つの渦糸を含むように，子午面上で渦面両側の平均子午線速度w m の方向に合わせるように行われる．下図は子午面上のメッシュを示したもので，点 O は修正済み，点m0 と点 P の位置 ∆m
m の間にはP は未修正，点 P0 は修正後の節点である．図を参照すれば点 P0 の位置 ∆m
次の関係が成立するすることが分かる．

m0
m
∆m
n ·∆m
=
wm
wm
n ·w
ただし n は i 一定線の単位法線ベクトルである．これより zik の修正値の式は次のようになる．
(n)

(n−1)

zik = zik(n−1)

0
+ α(zik
− zik

),

0
n ·∆m
m/n
n ·w
wm)
zik
= zi−1,k +wz (nただし α はダンピング係数である．rik の修正値の式も同様になる．なおこの計算は k 面が翼後縁で折れ曲らないように翼後縁から始められる．
line parallel toS »i-const line
»»»
»
S
¶
S
»
½S
S¶
»
½
»
k−1
S½ 0 S
S
1
P 
S
½
S
» wm
◦»
n
1


>
◦

:
»


½
S
» P
S½
»

k-surface in»A A mS»◦S∆mmeridian plane
k »O
A∆m
m0 S
S
¶
¶
»
S
»»S»
S
S »»
»
»
i
k+1
S
i−1
下流側の循環 Gmd = ΓD の値は Kutta の条件によって修正される．まず θ と tan βT の値が k 面の修正に伴って修正される．Kutta の条件は翼後縁における w = (wm , wu ) の方向が翼出口角に一致するときに満足されるから，ΓD の修正値の式は次のようになる．
(n)

(n−1)

ΓD = ΓD+ δΓD ,

δΓD = α2πr(wm tan βT −wu )44

非粘性流れの FEM 解析

周期境界 θ0 (z, r) は随伴渦面に一致するように，すなわち周期境界は各 k 面上で次式によって w に沿うよ
うに修正される．
(n)

θ0

(n−1)

= θ0

(n−1)

+ α(θ00 −θ0

),

1
wu
0
θ00 ≡ θi0k
= θi−1,0,k + ∆mrwm

k 面と周期境界の修正に伴って下流領域の要素分割が再度行われる．
下流境界における子午線速度 wmd = wmD は，一般化された釣り合い条件式 (11.63) を若干書換えた式
∆wm =

o
1 nω
1
2
∆ΓD −
∆Γ
D
2wm π
(2πr)2

によって修正される．更に下流境界の流量が上流境界の流量 GU になるように次式によって補正される．
(n)

(n−1)

wmD = wmD

− δwmD ,δwmD = P

´
³G
X
1
U
−
ρwmD r∆s
ρr∆s 2π

ΓD と wmD の修正に伴って下流境界の相対速度 w も修正される．
サブルーチン COEF： 有限要素方程式 (11.75) の係数行列すなわち剛性マトリックス

P

e
Hml
を設定する．

式 (11.74) の第 1 式に Gauss-Legendre 求積法を適用すればeHml=

Z 1 Z 1 Z 1³

8
3 ³
X
X
∂ψm ∂ψl ´∂ψm ∂ψl ´Jdξdη dζ =Jk∂xi ∂xi k−1 −1 −1 ∂xi ∂xi
i=1(11.84)k=1

ただし J はヤコビアン，添字 k は Gauss 点の値を意味する．剛性マトリックスは各要素ごと各 Gauss 点ごとに計算され重ね合わされる．始めに配列 nm に各要素の 8 つの節点番号が 11.3.2 項の図にしたがって入力される．要素剛性マトリックスは次の手順で計算される．式 (11.67) から psixi(1,l) = ψl,ξ , . . . の値，式 (11.69) からヤコビ行列 J の成分 aJm(i,j) の値，式 (11.68) から psix(i,l) = ∂ψl/∂xi とヤコビア
ン aJ = J の値が順に計算される．ただし sign は 11.3.2 項の表の値を記憶する配列，g は Gauss 点の座標である．また行列の反転にはサブルーチン MATINV が用いられる．ここではメモリー節約のため，有限要素方程式の係数行列 A はバンド幅の中だけが配列 H に記憶される．A の要素 amn は H の要素 H(m,l) に，特に対角要素 amm は H(m,0) に対応する．m = m，l = n−m．
次に境界条件について説明する．Dirichlet 条件の与えられる上流境界点では，Hml = 1 (l = 0), = 0 (l 6= 0)と置かれる．周期境界では，まず下方境界点の式が 相当の上方境界点の係数を重ね合わせることによって 完全なものに作り変えられる．それから上方境界点の式が 周期条件 (11.62c) の第 1 式になるように

Hml = 1 (l = 0), = −1 (l = −jf ), = 0 (その他の l) と置かれる．なお係数 Hml の添字 m はこの係
数が m 番目の式のものであること，また添字 l は φm+l の係数であることを示している．l の最小値はlb = −(kp +1)jp −1，その最大値は lf = (kp +1)jp +jf ，バンド幅は 2(kp +1)jp +jf +2 である．
サブルーチン CALPHI： 有限要素方程式 (11.75) の右辺を計算し，連立 1 次方程式を解いて速度ポテンシャ
ル φ を決定する．式 (11.74) の第 2 式は Gauss-Legendre 求積法を適用すれば
eKml=

Z 1Z 1Z 1
ψm ψl Jdξdη dζ =−1 −1 −1

8
X

(Jψm ψl )kk=1のように書換えられる．この右辺の値は各要素ごと各 Gauss 点ごとに入力される．Gauss 点 k の形状関数

psik = (ψl )k の値は式 (11.66) から計算され，J と b の値はそれぞれサブルーチン COEF または CALB で計算されたものが用いられる．なお質量の集中化 (mass lumping) はこのプログラムではしてはならない．

45

３次元翼列ポテンシャル流れのプログラム

右辺には Dirichlet 条件，Neumann 条件，周期条件が考慮される．Dirichlet 条件の課されている上流境
界点では，右辺は一般に 1−ΓU (θ−θ0 )/2π のように置くことができるが，ここでは ΓU = 0 ゆえ 1 になる．

Neumann 条件の課されている下流境界と翼表面では，右辺に −[絶対速度] · [面積要素ベクトル]，すなわち物理空間の境界要素から流入する流量が入力される．式 (11.74) の第 3 式は Gauss-Legendre 求積法を適用すれば次のようになる．
Z 1Z 14X
˜ dζ =
˜ m ψl )kGeml =ψm ψl Jdη
(Jψ
−1 −1

k=1

この式の ψm , ψl は立方体要素の形状関数としても良いが，ここでは境界における四辺形要素を正方形要素(−1 ≤ ξ, η ≤ 1) に写像しこの要素に対して定義された形状関数とする．この形状関数とその ξ, η に関する微分については 11.1.2 の式 (11.34) 以下を参照されたい．x 空間内の 4 節点 x 1 , x 2 , x 3 , x 4 の張る四辺形要素の面積ベクトル J̃ は次のようになる．
Ã! Ã ! x1xξ yξ zξψ ,ξ 
x2
J̃ ==xη yη zηψ ,η 
x3
x4

y1
y2
y3
y4


z1

z2 

z3 

z4

またこの要素から流入する流量は次のようになる．
¯¯¯ vx vy vz ¯¯¯¯¯
−vv · J̃ = − ¯ xξ yξ zξ ¯¯¯¯ xη yη zη ¯

(11.85)(11.86)下流境界の絶対速度 v の x, y, z 成分はその子午線成分 vm = wmD と周方向成分 vu = ΓD /(2πr) から求められ，また翼表面では n · w = 0，n · v = n · ω ×rr であるから，周速度 ω ×rr の y, z 成分が必要十分なものとなる．ハブとケーシング上では法線速度成分が 0 になるので何もする必要はなく，このことが FEM の利点のひとつにかぞえられている．周期境界では，下方境界点の右辺に上方境界点の右辺が加えられ，それから上方境界点の右辺がZ jF−
vu rdθ = −Γ/Nbladej0のように置かれる．
最後にこのようにして得られた連立 1 次方程式が Gauss 消去法で解かれ，速度ポテンシャル φ の値が決定される．この連立 1 次方程式の係数行列は COEF で入力されたものが用いられ，また Gauss 消去法の計算
は GAUSSB の中で行われる．
サブルーチン VELCTY： まず CALPHI で求めた φ を微分し相対速度 w を求め，それから相対 Mach 数 Mr ，
静圧 p などを求める．φ の微分は 11.3.5 項の方法で行われる．この方法で必要な 14 近傍節点は，図 11.5 に示すもので，あらかじめサブルーチン SETNEB で配列 neb に入力される．各節点ごとの局所近似 2 次多項式φ(x, y, z) は重み付き最小２乗法によって求められ，その計算は psi(l,k) = ψlk ，係数行列 a(n,l) = anl ，
その逆行列 a(l,n) = (a−1 )ln ，c(l) = cl の順に行われる．なお行列の反転はサブルーチン MATINV で行われる．絶対速度 v の xyz 方向成分は，vl = −cl , l = 1, 2, 3 となり，相対速度は次式から求められる．

wz = v1 ,wu = v2 cos θ+v3 sin θ−ωr,
wr = v2 sin θ−v3 cos θ46

非粘性流れの FEM 解析

また相対 Mach 数 Mr や静圧 p などの計算には次式が用いられる．

w2 (ωr)2
+
,
2
2
¡ ∗
(n)
(n−1)
(n−1) ¢∗
c2 = c2+ α c2 −c2,
c2 = (k−1)h,³ w2 ´1/2³ h ´k/(k−1)
c2
Mr =
,
h
=
,
p
=
p
,
0U
c2
k−1
h0U
h = I−

ρ=

k p
k−1 h

サブルーチン CALB： 式 (11.61) の右辺の b = b の値を計算する．比エンタルピー h の微分は，VELCTY における φ の微分と同様に行われる．∂h/∂xl = cl , l = 1, 2, 3．

b=

1³
∂h
∂h
∂h ´
w1
+w2
+w3
,
2
c
∂x1
∂x2
∂x3

w1 = wz ,w2 = wu cosθ+wr sinθ,

w3 = −wu sinθ+wr cosθ
プログラムの最後の部分は５つの補助的サブルーチンからなる．SETNEB は関数の局所近似に必要な 14 近
傍点の節点番号を配列 neb に入力するものである．固定壁境界における近傍点 k = 1, 2, . . . , 6 の配置換えは問題ないであろうか．残りの点 k = 7, 8, . . . , 14 の配置換えは隅の点を配慮して慎重に．なおここでは１
つ離れた節点に配置換えしている．WRTF2 と WRTF3 はそれぞれ子午面上の節点値または３次元計算領域のすべての節点値を出力するサブルーチンである．GAUSSB と MATINV は２次元翼列流れのものと同じである．
次に FORTRAN 95/90 Free source form で書かれた３次元翼列ポテンシャル流れのプログラムを示す．
PROGRAM MAIN!*****************************************************************************************! Problem: 3D Subsonic Potential Flow through a Mixed Flow Impeller! Numerical Methods: Variational FEM, Hexahedral Trilinear Elements, Gaussian Elimination!*****************************************************************************************DIMENSION z(0:20,0:6),th(0:20,0:6,0:6),r(0:20,0:6),x(3,1029),phi(1029), &w(3,1029),aM(1029),p(1029),tbt(0:6),Gmd(0:6),wmd(0:6),f(1029)
COMMON /COND/omg,pi,nob,Gmu,wmu,ak,aMu,rho0u /IJKP/na,i1,i2,k1,jp,kp,jkp &/DMAX/Drmax,Dthmax,DGmax,DMmax
CHARACTER*8 z1
jp=jf+1; kp=kf+1; jkp=jp*kp! ***** Input geometries of passage and vanes
DATA i1,i2,if,jf,k1,kf,mf/6,14,20,6,3,6,1029/
DATA z/-.300, -.160, -.020, .100, .200, .285, .350, .400, .465, .530, .600, &.670, .735, .800, .850, .890, .920, .950, .980, 1.010, 1.050,105*0., &-.278, -.230, -.180, -.130, -.080, -.040, .000, .050, .115, .180, .250, &.320, .385, .450, .500, .560, .630, .710, .800, .890, .963/
DATA r/ .250, .250, .250, .250, .252, .271, .300, .330, .369, .408, .450, &.492, .531, .570, .600, .622, .636, .648, .657, .665, .673,105*0., &.768, .763, .762, .767, .775, .786, .800, .820, .846, .872, .900, &.928, .954, .980, 1.000, 1.024, 1.052, 1.084, 1.119, 1.145, 1.160/
DATA th/6*0., 1.264, 1.169, 1.023, .879, .713, .528, .356, .167, 7*0.,105*0., &6*0., 1.264, 1.039, .814, .631, .457, .306, .181, .069, 7*0.,294*0., &6*0., 1.191, 1.088, .931, .782, .618, .453, .300, .138, 7*0.,105*0., &6*0., 1.191, .985, .788, .616, .451, .305, .183, .070, 7*0.,294*0., &6*0., 1.105, .996, .847, .700, .549, .396, .258, .117, 7*0.,105*0., &6*0., 1.105, .929, .756, .598, .448, .304, .186, .073, 7*0./! ***** Input flow conds
DATA omg,pi,nob,Gmu,wmu,ak,aMu,rho0u/1., 3.14159, 6, 0., .566, 1.4, .7, 1./

３次元翼列ポテンシャル流れのプログラム! ***** Discretize into elements, determine upstream bound conds and predict downstream!
bound conds
na=0
CALL SETX(z,th,r,x,tbt,Gmd,if,jf,kf,mf)
CALL PREDCT(z,th,r,w,wmd,if,jf,kf,mf);
GOTO 100! ***** Compute velocity potential phi, relative velocity w, Mach number aM and static!
pressure p
102 CALL CALB(th,r,x,w,if,jf,kf,mf)
100 na=na+1
CALL CORRCT(z,th,r,x,w,tbt,Gmd,wmd,if,jf,kf,mf)
CALL COEF(x,if,jf,kf,mf)
CALL CALPHI(z,th,r,x,phi, Gmd, wmd, if,jf,kf,mf)
CALL VELCTY(th,r,x,phi,w,aM,p,Gmd,if,jf,kf,mf)! ***** Output computational results
OPEN(20,FILE=’OUTPUT.DAT’)WRITE(20,60)na,Drmax,Dthmax,DGmax,DMmax
60 FORMAT(/1H ’na =’, I3, 2X ’Drmax =’, F7.4, 2X ’Dthmax=’, F7.4, &2X ’DGmax =’, F7.4, 2X ’DMmax =’, F7.4/)
IF(na<10)
GOTO 102
WRITE(20,61)tbt,Gmd,wmd
61 FORMAT(/1H 10X ’tbt =’, 7F7.3/1H0 10X ’Gmd =’, 7F7.3 &/1H0 10X ’wmd =’, 7F7.3/)
z1=’
z’; CALL WRTF2(z,z1,0)
z1=’r’; CALL WRTF2(r,z1,0)
DO k=1,7
SELECT CASE(k)CASE(1); z1=’ theta ’; FORALL(i=0:if,j=0:jf,k=0:kf)f(i*jkp+k*jp+j+1)=th(i,j,k)
CASE(2); z1=’
phi ’; FORALL(m=1:mf)f(m)=phi(m)CASE(3); z1=’
w_z ’; FORALL(m=1:mf)f(m)=w(1,m)
CASE(4); z1=’
w_u ’; FORALL(m=1:mf)f(m)=w(2,m)
CASE(5); z1=’
w_r ’; FORALL(m=1:mf)f(m)=w(3,m)
CASE(6); z1=’pressure’; FORALL(m=1:mf)f(m)=p(m)CASE(7); z1=’ Mach ’; FORALL(m=1:mf)f(m)=aM(m)ENDSELECTCALL WRTF3(f,z1,0)
ENDDO
CLOSE(20)! Drawing Mesh and static pressure distribution, Mach number contours, relative velocity! vectors on revolutional surfaces and meridian planes
CALL GRAPHICS(z,r,x,w,aM,p,if,jf,kf,mf)
STOP
END PROGRAM MAIN
! ********** Mesh Discretization
SUBROUTINE SETX(z,th,r,x,tbt,Gmd,if,jf,kf,mf)
DIMENSION z(0:if,0:kf),th(0:if,0:jf,0:kf),r(0:if,0:kf),x(3,mf),tbt(0:kf),Gmd(0:kf)
COMMON /COND/omg,pi,nob,Gmu,wmu,ak,aMu,rho0u /IJKP/na,i1,i2,k1,jp,kp,jkp &/SETX1/br(0:20),rt
DATA rt/.005/!rt=radius of trailing edge
pit=2.*pi/nob! ***** Generate mesh on meridional plane
DO 100 k=1,kf-1; t=k/FLOAT(kf); DO 100 i=0,ifz(i,k)=(1.-t)*z(i,0)+t*z(i,kf)r(i,k)=(1.-t)*r(i,0)+t*r(i,kf)
100 CONTINUE! ***** Compute breadths br of passage which are normal to k1-line
cycle_1: DO i=0,if

47

48

非粘性流れの FEM 解析

sz=z(i,kf)-z(i,0); sr=r(i,kf)-r(i,0)!span lengths
br(i)=SQRT(sz*sz+sr*sr)
IF(i==0.OR.i==if)CYCLE cycle_1
dz=z(i+1,k1)-z(i-1,k1); dr=r(i+1,k1)-r(i-1,k1)
dm=SQRT(dz*dz+dr*dr) ; zn=-dr/dm; rn=dz/dm!unit normal of dm
br(i)=zn*sz+rn*sr!passage breadths
ENDDO cycle_1! ***** Generate mesh on blade surfaces
DO k=1,kf-1; t=k/FLOAT(kf)c0=(1.-t)*(1.-2.*t); c1=4.*t*(1.-t); c2=-t*(1.-2.*t)
DO i=i1,i2th(i, 0,k)=c0*th(i, 0,0)+c1*th(i, 0,k1)+c2*th(i, 0,kf)th(i,jf,k)=c0*th(i,jf,0)+c1*th(i,jf,k1)+c2*th(i,jf,kf)
ENDDO
ENDDO
! ***** Predict downstream bound conds GMD
i=i2
DO k=0,kf
nn=0
101 nn=nn+1; cbt=1./SQRT(1.+tbt(k)*tbt(k)); sbt=cbt*tbt(k)dz=z(i,k)-z(i-1,k); rm=(r(i-1,k)+r(i,k))/2.
dr=r(i,k)-r(i-1,k); dm=SQRT(dz*dz+dr*dr)
dth=th(i, 0,k)-th(i-1, 0,k); tbt0=(rm*dth+rt*cbt)/(dm-rt*sbt); bt0=ATAN(tbt0)
dth=th(i,jf,k)-TH(i-1,jf,k); tbt1=(rm*dth-rt*cbt)/(dm+rt*sbt); bt1=ATAN(tbt1)
tbt(k)=TAN((bt0+bt1)/2.); IF(nn<3)
GOTO 101
rr=(r(i,k)-r(i,k1))*2./(r(i,kf)-r(i,0))
wm=wmu*r(0,k1)*br(0)/r(i,k1)/br(i)*(1.+.2*rr)!suppose semi-forced vortex design
Gmd(k)=2.*pi*r(i,k)*(wm*tbt(k)+omg*r(i,k))
ENDDO! ***** Generate mesh on periodic boundaries
DO k=0,kf
DO i=i1-1,0,-1
dz=z(i+1,k)-z(i,k); rm=(r(i,k)+r(i+1,k))/2.
dr=r(i+1,k)-r(i,k); dm=SQRT(dz*dz+dr*dr)
wu=Gmu/(2.*pi*rm)-omg*rm
wm=wmu*r(0,k1)*br(0)/r(i,k1)/br(i)th(i,0,k)=th(i+1,0,k)-wu/wm*dm/rm
ENDDO
DO i=i2+1,if
dz=z(i,k)-z(i-1,k); rm=(r(i-1,k)+r(i,k))/2.
dr=r(i,k)-r(i-1,k); dm=SQRT(dz*dz+dr*dr)
wu=Gmd(k)/(2.*pi*rm)-omg*rm
wm=wmu*r(0,k1)*br(0)/r(i,k1)/br(i)th(i,0,k)=th(i-1,0,k)+wu/wm*dm/rm
ENDDO
ENDDO
! ***** Determine coordinates of nodal points
DO 102 i=0,if; DO 102 k=0,kf
IF(i<i1.OR.i>i2)th(i,jf,k)=th(i,0,k)th(i,jf,k)=th(i,jf,k)+pit
DO j=0,jf; t=j/FLOAT(jf); m=i*jkp+k*jp+j+1
th(i,j,k)=(1.-t)*th(i,0,k)+t*th(i,jf,k)x(1,m)=z(i,k)
x(2,m)=r(i,k)*SIN(th(i,j,k))x(3,m)=r(i,k)*COS(th(i,j,k))
ENDDO
102 CONTINUE

３次元翼列ポテンシャル流れのプログラム
END SUBROUTINE SETX
! ********** Determine Upstream Bound Conds and Predict Downstream Bound Conds
SUBROUTINE PREDCT(z,th,r,w,wmd,if,jf,kf,mf)
DIMENSION z(0:if,0:kf),th(0:if,0:jf,0:kf),r(0:if,0:kf),w(3,mf),wmd(0:kf)
COMMON /COND/omg,pi,nob,Gmu,wmu,ak,aMu,rho0u /IJKP/na,i1,i2,k1,jp,kp,jkp &/SETX1/br(0:20),rt /PRED1/rhou,h0u,rot,p0u,flw! ***** Determine upstream bound conds
vuu=Gmu/(2.*pi*r(0,k1));
wuu=vuu-omg*r(0,k1)
wu =SQRT(wmu*wmu+wuu*wuu); cu =wu/aMu
hu =cu*cu/(ak-1.);h0u=hu+(wmu*wmu+vuu*vuu)/2.
!stagnation enthalpy
rot=hu+wu*wu/2.-omg*omg*r(0,k1)*r(0,k1)/2.
!rothalpyp0u=(ak-1.)/ak*rho0u*h0u!stagnation pressure
flw=0.
DO k=0,kf-1
DZ=Z(0,K+1)-Z(0,K); RM=(R(0,K)+R(0,K+1))/2.
DR=R(0,K+1)-R(0,K); DS=SQRT(DZ*DZ+DR*DR)VUU=GMU/(2.*PI*RM); HU=H0U-(WMU*WMU+VUU*VUU)/2.
RHO1=RHO0U*(HU/H0U)**(1./(AK-1.))
IF(K.EQ.K1)RHOU=RHO1!density at mid-passage
FLW=FLW+2.*PI*RHO1*WMU*RM*DS!mass flow rate
ENDDO
DO 110 k=0,kf; DO 110 j=0,jf; m=k*jp+j+1w(1,m)=wmu*(r(0,kf)-r(0,0))/br(0)w(2,m)=Gmu/(2.*pi*r(0,k))-omg*r(0,k)
w(3,m)=-wmu*(z(0,kf)-z(0,0))/br(0)110 CONTINUE! ***** Predict downstream bound conds
i=if; Swmd=0.
cycle_2: DO k=0,kf
rr=(r(i,k)-r(i,k1))*2./(r(i,kf)-r(i,0))
wmd(k)=wmu*r(0,k1)*br(0)/r(i,k1)/br(i)*(1.+.2*rr)!semi-forced vortex
IF(k==0)CYCLE cycle_2
dz=z(i,k)-z(i,k-1); rm=(r(i,k-1)+r(i,k))/2.
dr=r(i,k)-r(i,k-1); ds=SQRT(dz*dz+dr*dr)
Swmd=Swmd+(wmd(k-1)+wmd(k))/2.*rm*ds
ENDDO cycle_2
dwmd=(flw/(2.*pi*rhou)-Swmd)/(r(i,k1)*br(i))!flow compensation
FORALL(k=0:kf)wmd(k)=wmd(k)+dwmd
END SUBROUTINE PREDCT! ********** Correct Downstream Region
SUBROUTINE CORRCT(z,th,r,x,w,tbt,Gmd,wmd,if,jf,kf,mf)
DIMENSION z(0:if,0:kf),th(0:if,0:jf,0:kf),r(0:if,0:kf),x(3,mf),w(3,mf), &tbt(0:kf),Gmd(0:kf),wmd(0:kf),w1(0:6),w2(0:6),w3(0:6),w4(0:6),al3(20)COMMON /COND/omg,pi,nob,Gmu,wmu,ak,aMu,rho0u /IJKP/na,i1,i2,k1,jp,kp,jkp &/DMAX/Drmax,Dthmax,DGmax,DMmax /SETX1/br(0:20),rt &/PRED1/rhou,h0u, rot,p0u,flw /VELC2/rho(1029)
DATA al1,al2,al3/.6, .2, 2*.25, .5, 17*1./
Drmax=0.; Dthmax=0.; DGmax=0.; pit=2.*pi/nob
IF(na==1)
GOTO 200! ***** Correct k-surface in order to coincide with a vortex filament
DO k=1,kf-1; DO i=i2,if; m=i*jkp+k*jp+1
wz=w(1,m)+w(1,m+jf)+w(1,m-jkp)+w(1,m-jkp+jf)
wr=w(3,m)+w(3,m+jf)+w(3,m-jkp)+w(3,m-jkp+jf)
dz=z(i,k)-z(i-1,k); rn= z(i,kf)-z(i,0)
dr=r(i,k)-r(i-1,k); zn=-r(i,kf)+r(i,0)!normal of span

49

50

非粘性流れの FEM 解析

wn=zn*wz+rn*wr; dmn=zn*dz+rn*dr!normal componentsz1=z(i-1,k)+wz*dmn/wn; Dz=al1*(z1-z(i,k))!correction pts on span
r1=r(i-1,k)+wr*dmn/wn; Dr=al1*(r1-r(i,k))
Drmax=AMAX1(Drmax,ABS(Dr))
z(i,k)=z(i,k)+Dzr(i,k)=r(i,k)+Dr
ENDDO; ENDDO! ***** Correct circulation Gamma_D by the Katta condition
i=i2; j=jf
FORALL(k=0:kf)
w1(k)=z(i,k); w2(k)=th(i,0,k); w3(k)=r(i,k); w4(k)=tbt(k)ENDFORALL
DO k=1,kf-1; dr=w3(k+1)-w3(k)!quadr interp along trailing edge
al=(r(i,k)-w3(k))/dr; am=(w3(k)-w3(k-1))/dr
c0=al*(al-1.)/am/(am+1.)
c1=-(al+am)*(al-1.)/am
c2=(al+am)*al/(am+1.)th(i,0,k)=c0*w2(k-1)+c1*w2(k)+c2*w2(k+1)
tbt(k)=c0*w4(k-1)+c1*w4(k)+c2*w4(k+1)
ENDDO
DO k=0,kf; m=i*jkp+k*jp+1
wz=(w(1,m)+w(1,m+jf))/2.; wu=(w(2,m)+w(2,m+jf))/2.
wr=(w(3,m)+w(3,m+jf))/2.; wm=SQRT(wz*wz+wr*wr)
DGmd=al2*2.*pi*r(i,k)*(wm*tbt(k)-wu)
DGmax=AMAX1(DGmax,ABS(DGmd))
Gmd(k)=Gmd(k)+DGmd
ENDDO! ***** correct trailing vortex sheet in order to coincide with vortex filaments
DO i=i2+1,if; DO k=0,kf; m=i*jkp+k*jp+1
wz=w(1,m)+w(1,m+jf)+w(1,m-jkp)+w(1,m-jkp+jf)
wu=w(2,m)+w(2,m+jf)+w(2,m-jkp)+w(2,m-jkp+jf)
wr=w(3,m)+w(3,m+jf)+w(3,m-jkp)+w(3,m-jkp+jf)
dz=z(i,k)-z(i-1,k); rm=(r(i-1,k)+r(i,k))/2.
dr=r(i,k)-r(i-1,k); dm=SQRT(dz*dz+dr*dr)
wm=SQRT(wz*wz+wr*wr)
th1=th(i-1,0,k)+wu/wm*dm/rm; Dth=al3(na)*(th1-th(i,0,k))
Dthmax=AMAX1(Dthmax,ABS(Dth))th(i,0,k)=th(i,0,k)+Dth
ENDDO; ENDDO
DO 213 i=i2,if; DO 213 j=0,jf; DO 213 k=0,kf; m=i*jkp+k*jp+j+1
th(i,j,k)=th(i,0,k)+(j/FLOAT(jf))*pitx(1,m)=z(i,k)
x(2,m)=r(i,k)*SIN(th(i,j,k))x(3,m)=r(i,k)*COS(th(i,j,k))
213 CONTINUE! ***** Correct downstream meridian velocities w_mD
200 nn=0; nnmax=10; i=if
220 nn=nn+1; Srhowm=0.; Srho=0.
DO k=k1-1,0,-1
rm=(r(i,k)+r(i,k+1))/2.; g0=Gmd(k); g1=Gmd(k+1)
wmd(k)=wmd(k+1)-(omg/pi*(g1-g0)-(g1*g1-g0*g0)/(4.*pi*pi*rm*rm))/(wmd(k)+wmd(k+1))
ENDDO
DO k=k1+1,kf
rm=(r(i,k-1)+r(i,k))/2.; g0=Gmd(k); g1=Gmd(k-1)
wmd(k)=wmd(k-1)+(omg/pi*(g0-g1)-(g0*g0-g1*g1)/(4.*pi*pi*rm*rm))/(wmd(k-1)+wmd(k))
ENDDO
DO k=1,kf; m=if*jkp+(k+1)*jp

３次元翼列ポテンシャル流れのプログラム
dz=z(i,k)-z(i,k-1); rm=(r(i,k-1)+r(i,k))/2.
dr=r(i,k)-r(i,k-1); ds=SQRT(dz*dz+dr*dr)
rhom=(rho(m-jp)+rho(m))/2.; IF(na.EQ.1)rhom=rhou
Srhowm=Srhowm+rhom*(wmd(k-1)+wmd(k))/2.*rm*ds
Srho =Srho +rhom*rm*ds
ENDDO
Dwmd=(flw/(2.*pi)-Srhowm)/Srho
FORALL(k=0:kf)wmd(k)=wmd(k)+Dwmd
IF(nn>nnmax)STOP ’difficult to obtain meridian velocity at downstream boundary’IF(ABS(Dwmd)>.001.OR.nn<3)
GOTO 220
DO k=0,kf; DO j=0,jf; m=if*jkp+k*jp+j+1w(1,m)= wmd(k)*(r(i,kf)-r(i,0))/br(i)w(2,m)= Gmd(k)/(2.*pi*r(i,k))-omg*r(i,k)
w(3,m)=-wmd(k)*(z(i,kf)-z(i,0))/br(i)ENDDO; ENDDO
END SUBROUTINE CORRCT! ********** Set up Coefficient Matrix of Linear Equations
SUBROUTINE COEF(x,if,jf,kf,mf)
DIMENSION x(3,mf),psixi(3,8),aJm(3,3),psix(3,8),sign(3,8),aw(3),iw(3)COMMON /IJKP/na,i1,i2,k1,jp,kp,jkp /COEF1/H(1029,-57:62),aJ(8,720), &nm(8,720),lb,lf,nf
lb=-jkp-jp-1; lf=jkp+jp+jf; n0=jf*kf; nf=if*jf*kf! ***** Determine node numbers of each elements
n=0; m=0; 300 n=n+1; m=m+1
nm(1,n)=m;nm(2,n)=m+jkp;
nm(3,n)=m+1;
nm(4,n)=m+jp
nm(5,n)=m+jp+1; nm(6,n)=m+jkp+jp; nm(7,n)=m+jkp+1; nm(8,n)=nm(6,n)+1IF(MOD(n,jf)==0)m=m+1
IF(MOD(n,n0)==0)m=m+jp; IF(n<nf)
GOTO 300! ***** Set up coefficient matrix H
g=1./1.732
DATA sign/-1,-1,-1, 1,-1,-1, -1, 1,-1, -1,-1, 1, &-1, 1, 1, 1,-1, 1, 1, 1,-1, 1, 1, 1/
FORALL(m=1:mf,l=lb:lf)H(m,l)=0.
DO n=1,nf; DO k=1,8
DO l=1,8
psixi(1,l)=sign(1,l)*(1.+sign(2,l)*sign(2,k)*g)*(1.+sign(3,l)*sign(3,k)*g)/8.
psixi(2,l)=sign(2,l)*(1.+sign(3,l)*sign(3,k)*g)*(1.+sign(1,l)*sign(1,k)*g)/8.
psixi(3,l)=sign(3,l)*(1.+sign(1,l)*sign(1,k)*g)*(1.+sign(2,l)*sign(2,k)*g)/8.
ENDDO
DO i=1,3; DO j=1,3; aJm(i,j)=0.; DO l=1,8aJm(i,j)=aJm(i,j)+psixi(i,l)*x(j,nm(l,n))
ENDDO; ENDDO; ENDDO
aJ(k,n)=1.; e=1.E-5
CALL MATINV(aJm,aw,iw,3,aJ(k,n),e)DO i=1,3; DO l=1,8; psix(i,l)=0.; DO j=1,3
psix(i,l)=psix(i,l)+aJm(i,j)*psixi(j,l)
ENDDO; ENDDO; ENDDO
DO me=1,8; m=nm(me,n)
DO le=1,8; l=nm(le,n)-m; DO i=1,3
H(m,l)=H(m,l)+aJ(k,n)*psix(i,me)*psix(i,le)
ENDDO; ENDDO; ENDDO
ENDDO; ENDDO! ***** Take account of Dirichlet bound cond
FORALL(m=1:jkp,l=lb:lf)H(m,l)=0.
!upstream bound
FORALL(m=1:jkp)H(m,0)=1.
! ***** Take account of periodic bound conds

51

52

非粘性流れの FEM 解析

cycle_3: DO i=1,if; IF(i>i1.AND.i<i2)CYCLE cycle_3
DO k=0,kf; m=i*jkp+k*jp+1
DO 303 ii=-1,1; DO 303 jj=-1,0; DO 303 kk=-1,1IF(i+ii<0.OR.i+ii>if)
GOTO 303
IF(k+kk<0.OR.k+kk>kf)
GOTO 303
l=ii*jkp+kk*jp+jj
H(m,l+jf)=H(m+jf,l)
303 CONTINUEENDDOENDDO cycle_3
cycle_4: DO i=1,if; IF(i>i1.AND.i<i2)CYCLE cycle_4
DO k=0,kf; m=i*jkp+(k+1)*jp
FORALL(l=lb:lf)H(m,l)=0.
H(m,0)=1.; H(m,-jf)=-1.
ENDDO
ENDDO cycle_4
END SUBROUTINE COEF

!periodic bound j=0!periodic bound jf! ********** Solve Finte Element Eqns of phi by Gaussian Elimination
SUBROUTINE CALPHI(z,th,r,x,phi,Gmd,wmd,if,jf,kf,mf)
DIMENSION z(0:if,0:kf),TH(0:if,0:jf,0:kf),r(0:if,0:kf),x(3,mf),phi(mf),Gmd(0:kf),wmd(0:kf), &sign(3,8),psik(8),nmd(4,36),nmb(4,48),v(3,0:6),signb(2,4),psikb(4),psixi(2,4),tJ(2,3)
COMMON /COND/omg,pi,nob,Gmu,wmu,ak,aMu,rho0u /IJKP/na,i1,i2,k1,jp,kp,jkp &/COEF1/H(1029,-57:62),aJ(8,720),nm(8,720),lb,lf,nf /CALB1/b(1029)
n0=jf*kf; nb=(i2-i1)*kf! ***** Set up rhs of linear eqns
g=1./1.732
DATA sign/-1,-1,-1, 1,-1,-1, -1, 1,-1, -1,-1, 1, &-1, 1, 1, 1,-1, 1, 1, 1,-1, 1, 1, 1/
FORALL(m=1:mf)phi(m)=0.
DO n=1,nf; DO k=1,8
FORALL(l=1:8)psik(l)=(1.+sign(1,l)*sign(1,k)*g)*(1.+sign(2,l)*sign(2,k)*g) &*(1.+sign(3,l)*sign(3,k)*g)/8.
DO me=1,8; m=nm(me,n); DO le=1,8; l=nm(le,n)phi(m)=phi(m)-aJ(k,n)*psik(me)*psik(le)*b(l)ENDDO; ENDDO
ENDDO; ENDDO! ***** Take account of Dirichlet condition
DO k=0,kf; DO j=0,jf; m=k*jp+j+1!upstream bc
phi(m)=1.-Gmu/(2.*pi)*(th(0,j,k)-th(0,0,0))
ENDDO; ENDDO! ***** Take account of Neumann condition at downstream boundary
n=0; m=if*jkp; 310 n=n+1; m=m+1
nmd(1,n)=m; nmd(2,n)=m+1; nmd(3,n)=m+jp; nmd(4,n)=m+jp+1
IF(MOD(n,jf)==0)m=m+1; IF(n<n0)
GOTO 310
DO k=0,kf
dz=z(if,k)-z(if-1,k); dr=r(if,k)-r(if-1,k); dm=SQRT(dz*dz+dr*dr)v(1,k)=(dz/dm)*wmd(k)v(2,k)=Gmd(k)/(2.*pi*r(if,k))v(3,k)=(dr/dm)*wmd(k)ENDDO
DATA signb/-1,-1, 1,-1, -1,1, 1,1/
DO n=1,n0; DO k=1,4
FORALL(l=1:4)
psikb(l)=(1.+signb(1,l)*signb(1,k)*g)*(1.+signb(2,l)*signb(2,k)*g)/4.
psixi(1,l)=signb(1,l)*(1.+signb(2,l)*signb(2,k)*g)/4.
psixi(2,l)=signb(2,l)*(1.+signb(1,l)*signb(1,k)*g)/4.

３次元翼列ポテンシャル流れのプログラム
ENDFORALLFORALL(i=1:2,j=1:3)tJ(i,j)=0.
DO j=1,3; DO l=1,4tJ(1,j)=tJ(1,j)+psixi(1,l)*x(j,nmd(l,n))tJ(2,j)=tJ(2,j)+psixi(2,l)*x(j,nmd(l,n))
ENDDO; ENDDO
DO le=1,4; l=nmd(le,n); j0=MOD(l-1,jp); k0=MOD((l-1)/jp,kp)
costh=COS(th(if,j0,k0)); sinth=SIN(th(if,j0,k0))v1=v(1,k0); v2= v(2,k0)*costh+v(3,k0)*sinth
v3=-v(2,k0)*sinth+v(3,k0)*costh
vtJ=v1*(tJ(1,2)*tJ(2,3)-tJ(2,2)*tJ(1,3)) &+v2*(tJ(1,3)*tJ(2,1)-tJ(2,3)*tJ(1,1)) &+v3*(tJ(1,1)*tJ(2,2)-tJ(2,1)*tJ(1,2))
DO me=1,4; m=nmd(me,n)phi(m)=phi(m)-psikb(me)*psikb(le)*vtJ
ENDDO
ENDDO
ENDDO; ENDDO! ***** Take account of Neumann condition on blade surfaces
n=0; m=i1*jkp+1-jp; 311 n=n+1; m=m+jp
nmb(1,n)=m; nmb(2,n)=m+jkp; nmb(3,n)=m+jp; nmb(4,n)=m+jkp+jp
IF(MOD(n,kf)==0)m=m+jp; IF(n<nb)
GOTO 311
j0=0; sign0=-1.
320 CONTINUE
DO n=1,nb; DO k=1,4
FORALL(l=1:4)
psikb(l)=(1.+signb(1,l)*signb(1,k)*g)*(1.+signb(2,l)*signb(2,k)*g)/4.
psixi(1,l)=signb(1,l)*(1.+signb(2,l)*signb(2,k)*g)/4.
psixi(2,l)=signb(2,l)*(1.+signb(1,l)*signb(1,k)*g)/4.
ENDFORALLFORALL(i=1:2,j=1:3)tJ(i,j)=0.
DO j=1,3; DO l=1,4tJ(1,j)=tJ(1,j)+psixi(1,l)*x(j,nmb(l,n))tJ(2,j)=tJ(2,j)+psixi(2,l)*x(j,nmb(l,n))
ENDDO; ENDDO
DO le=1,4; l=nmb(le,n); i0=(l-1)/jkp; k0=MOD((l-1)/jp,kp)
v2= omg*r(i0,k0)*COS(th(i0,j0,k0))v3=-omg*r(i0,k0)*SIN(th(i0,j0,k0))
vtJ=v2*(tJ(1,3)*tJ(2,1)-tJ(2,3)*tJ(1,1))+v3*(tJ(1,1)*tJ(2,2)-tJ(2,1)*tJ(1,2))
DO me=1,4; m=nmb(me,n)phi(m)=phi(m)+sign0*psikb(me)*psikb(le)*vtJ
ENDDO; ENDDO
ENDDO; ENDDOFORALL(l=1:4,n=1:nb)nmb(l,n)=nmb(l,n)+jf
j0=j0+jf; sign0=1.; IF(j0==jf)
GOTO 320! ***** Take account of periodic conds
cycle_5: DO i=1,if; IF(i>i1.AND.i<i2)CYCLE cycle_5
DO k=0,kf; m=i*jkp+k*jp+1
phi(m)=phi(m)+phi(m+jf)!periodic bc j=0
m=m+jf; IF(i<=i1)phi(m)=-Gmu/nob!periodic bc jf
IF(i>=i2)phi(m)=-Gmd(k)/nob!periodic bc jf
ENDDO
ENDDO cycle_5! ***** Determine potential phi solving linear eqns by Gaussian elimination
CALL GAUSSB(H,phi,mf,lb,lf,det)
END SUBROUTINE CALPHI

53

54

非粘性流れの FEM 解析! ********** Compute Velocities
SUBROUTINE VELCTY(th,r,x,phi,w,aM,p,Gmd,if,jf,kf,mf)
DIMENSION th(0:if,0:jf,0:kf),r(0:if,0:kf),x(3,mf),phi(mf),w(3,mf),aM(mf),p(mf), &Gmd(0:kf),neb(14,1029),bphi(14),psi(9,14),c(3),a(9,9),aw(9),iw(9)COMMON /COND/omg,pi,nob,Gmu,wmu,ak,aMu,who0u /IJKP/na,i1,i2,k1,jp,kp,jkp &/DMAX/Drmax,Dthmax,DGmax,DMmax /VELC1/cc(1029),h(1029) /VELC2/rho(1029) &/SETX1/br(0:20),rt /PRED1/rhou,h0u,rot,p0u,flw
DATA al5/.25/
DMmax=0.; pit=2.*pi/nob! ***** Determine neighbouring point numbers of each nodes
CALL SETNEB(neb,if,jf,kf,mf)
OPEN(20,FILE=’OUTPUT.DAT’)! ***** Compute relative velocities w
DO m=1+jkp,mf-jkp; m0=m!0<i<if
i=(m-1)/jkp; j=MOD(m-1,jp); k0=MOD((m-1)/jp,kp)
DO k=1,14; mk=neb(k,m0)
ki=(mk-1)/jkp; kj=MOD(mk-1,jp); kk=MOD((mk-1)/jp,kp)
bphi(k)=phi(mk)-phi(m0)
FORALL(l=1:3)psi(l,k)=x(l,mk)-x(l,m0)
IF(i>=i1.AND.i<=i2)THENELSEIF(j==0.AND.kj==jf-1)THEN
IF(i<i1)bphi(k)=bphi(k)+Gmu/nob!periodic bound j=0
IF(i>i2)bphi(k)=bphi(k)+Gmd(kk)/nobpsi(2,k)=r(ki,kk)*SIN(th(ki,kj,kk)-pit)-x(2,m0)
psi(3,k)=r(ki,kk)*COS(th(ki,kj,kk)-pit)-x(3,m0)
ELSEIF(j==jf.AND.kj==1)THEN
IF(i<i1)bphi(k)=bphi(k)-Gmu/nob!periodic bound jf
IF(i>i2)bphi(k)=bphi(k)-Gmd(kk)/nobpsi(2,k)=r(ki,kk)*SIN(th(ki,kj,kk)+pit)-x(2,m0)
psi(3,k)=r(ki,kk)*COS(th(ki,kj,kk)+pit)-x(3,m0)
ENDIF
DO l1=1,3; l2=MOD(l1,3)+1; l3=MOD(l2,3)+1
psi(l1+3,k)=psi(l1,k)*psi(l1,k)psi(l1+6,k)=2.*psi(l2,k)*psi(l3,k)
ENDDO
ENDDO
DO 401 n=1,9; DO 401 l=1,9; a(n,l)=0.; DO 401 k=1,14
wt=psi(4,k)+psi(5,k)+psi(6,k); wt=1./wt/wt/wt
a(n,l)=a(n,l)+wt*psi(n,k)*psi(l,k)
401 CONTINUE
nf=9; e=1.E-5
CALL MATINV(a,aw,iw,nf,det,e)
DO 402 l=1,3; c(l)=0.; DO 402 n=1,9; DO 402 k=1,14
wt=psi(4,k)+psi(5,k)+psi(6,k); wt=1./wt/wt/wt
c(l)=c(l)+a(l,n)*wt*psi(n,k)*bphi(k)402 CONTINUE
costh=COS(th(i,j,k0)); sinth=SIN(th(i,j,k0))w(1,m)=-c(1)!w_zw(2,m)=-c(2)*costh+c(3)*sinth-omg*r(i,k0)!w_u
w(3,m)=-c(2)*sinth-c(3)*costh
!w_r
ENDDO
! ***** Compute Mach number aM and static pressure p
DO m=1,mf; i=(m-1)/jkp; k0=MOD((m-1)/jp,kp)
ww=w(1,m)*w(1,m)+w(2,m)*w(2,m)+w(3,m)*w(3,m)
cc1=(ak-1.)*(rot-ww/2.+omg*omg*r(i,k0)*r(i,k0)/2.)IF(na==1)THEN; cc(m)=cc1
ELSE; cc(m) =cc(m)+al5*(cc1-cc(m)); ENDIF

３次元翼列ポテンシャル流れのプログラム
aM1=aM(m)
aM(m) =SQRT(ww/cc(m))!relative Mach number
DMmax =AMAX1(DMmax,ABS(aM1-aM(m)))h(m) =cc(m)/(ak-1.); ak1=ak/(ak-1.)!static enthalpy
p(m) =p0u*(h(m)/h0u)**ak1!static pressurerho(m)=ak1*p(m)/h(m)!densityIF(na==1)rho(m)=rhou+.5*(rho(m)-rhou)
ENDDO! ***** Take mean values of p at leading and trailing edges
DO k=0,kf; m=i1*jkp+k*jp+1
p(m)=(p(m)+p(m+jf))/2.; p(m+jf)=p(m)m=i2*jkp+k*jp+1
p(m)=(p(m)+p(m+jf))/2.; p(m+jf)=p(m)ENDDO
END SUBROUTINE VELCTY! ********** Compute rhs of Linear Eqns
SUBROUTINE CALB(th,r,x,w,if,jf,kf,mf)
DIMENSION th(0:if,0:jf,0:kf),r(0:if,0:kf),x(3,mf),w(3,mf),neb(14,1029),bh(14), &psi(9,14),c(3),a(9,9),aw(9),iw(9)COMMON /COND/omg,pi,nob,Gmu,wmu,ak,aMu,rho0u /IJKP/na,i1,i2,k1,jp,kp,jkp &/CALB1/b(1029) /VELC1/cc(1029),h(1029)
CALL SETNEB(neb,if,jf,kf,mf)
pit=2.*pi/nob
DO m=1+jkp,mf-jkp; m0=mi=(m-1)/jkp; j=MOD(m-1,jp); k0=MOD((m-1)/jp,kp)
DO k=1,14; mk=neb(k,m0)
ki=(mk-1)/jkp; kj=MOD(mk-1,jp); kk=MOD((mk-1)/jp,kp)
bh(k)=h(mk)-h(m0)
FORALL(l=1:3)psi(l,k)=x(l,mk)-x(l,m0)
IF(i>=i1.AND.i<=i2)THENELSEIF(j==0.AND.kj==jf-1)THEN
psi(2,k)=r(ki,kk)*SIN(th(ki,kj,kk)-pit)-x(2,m0)!periodic bound j=0
psi(3,k)=r(ki,kk)*COS(th(ki,kj,kk)-pit)-x(3,m0)
ELSEIF(j==jf.AND.kj==1)THEN
psi(2,k)=r(ki,kk)*SIN(th(ki,kj,kk)+pit)-x(2,m0)!periodic bound jf
psi(3,k)=r(ki,kk)*COS(th(ki,kj,kk)+pit)-x(3,m0)
ENDIF
DO l1=1,3; l2=MOD(l1,3)+1; l3=MOD(l2,3)+1
psi(l1+3,k)=psi(l1,k)*psi(l1,k)psi(l1+6,k)=2.*psi(l2,k)*psi(l3,k)
ENDDO
ENDDO
DO 411 n=1,9; DO 411 l=1,9; a(n,l)=0.; DO 411 k=1,14
wt=psi(4,k)+psi(5,k)+psi(6,k); wt=1./wt/wt/wt
a(n,l)=a(n,l)+wt*psi(n,k)*psi(l,k)
411 CONTINUE
nf=9; e=1.E-5
CALL MATINV(a,aw,iw,nf,det,e)
DO 412 l=1,3; c(l)=0.; DO 412 n=1,9; DO 412 k=1,14
wt=psi(4,k)+psi(5,k)+psi(6,k); wt=1./wt/wt/wt
c(l)=c(l)+a(l,n)*wt*psi(n,k)*bh(k)412 CONTINUE
costh=COS(th(i,j,k0)); sinth=SIN(th(i,j,k0))
w2= w(2,m)*costh+w(3,m)*sinth
w3=-w(2,m)*sinth+w(3,m)*costhb(m)=(w(1,m)*c(1)+w2*c(2)+w3*c(3))/cc(m)

55

56

非粘性流れの FEM 解析

ENDDO
END SUBROUTINE CALB
! ********** SUBSIDIARY SUBROUTINES! ********** Determine Neighbouring Nodal Points
SUBROUTINE SETNEB(neb,if,jf,kf,mf)
INTEGER neb(14,mf)
COMMON /IJKP/na,i1,i2,k1,jp,kp,jkp
DO 10 k=1,14
SELECT CASE(k)CASE(1,7,8,11,12); ii=-1; CASE(3:6); ii=0; CASE DEFAULT; ii=1
ENDSELECTSELECT CASE(k)CASE(3,7:10); ij=-1; CASE(1,2,5,6); ij=0; CASE DEFAULT; ij=1
ENDSELECTSELECT CASE(k)CASE(5,7,9,11,13); ik=-1; CASE(:4); ik=0; CASE DEFAULT; ik=1
ENDSELECT
DO m=1+jkp,mf-jkp!0<i<ifneb(k,m)=m+ii*jkp+ij+ik*jp
ENDDO
10 CONTINUE
DO 11 i=1,if-1; DO 11 j=0,jf; m=i*jkp+j+1
neb( 5,m)=m+3*jp;
neb( 7,m)=neb(5,m)-1!on hub
neb( 9,m)=neb(5,m)-jkp; neb(11,m)=neb(5,m)+1
neb(13,m)=neb(5,m)+jkp; m=m+jp*kf
neb( 6,m)=m-3*jp;
neb( 8,m)=neb(6,m)-1!on casingneb(10,m)=neb(6,m)-jkp; neb(12,m)=neb(6,m)+1
neb(14,m)=neb(6,m)+jkp
11 CONTINUE
DO 12 i=i1,i2; DO 12 k0=0,kf; m=i*jkp+k0*jp+1!on bladeneb(3,m)=m+3; neb( 7,m)=neb(3,M)-jkp!pressure side
neb( 8,m)=neb(3,M)+jkp
IF(k0/= 0)neb( 9,m)=neb(3,M)-jp
IF(k0/=kf)
neb(10,m)=neb(3,M)+jp; m=m+jf
neb(4,m)=m-3; neb(11,m)=neb(4,M)-jkp!suction sideneb(12,m)=neb(4,M)+jkp
IF(k0/= 0)neb(13,m)=neb(4,M)-jp
IF(k0/=kf)
neb(14,m)=neb(4,M)+jp
12 CONTINUE
cycle_1: DO i=1,if-1!periodic boundIF(i>=i1.AND.i<=i2)CYCLE cycle_1
DO k0=0,kf; DO k=1,14; m=i*jkp+k0*jp+1
IF(MOD(neb(k,m),jp)==0)neb(k,m)=neb(k,m)+jf; m=m+jf!j=0
IF(MOD(neb(k,m),jp)==1)neb(k,m)=neb(k,m)-jf!jf
ENDDO; ENDDO
ENDDO cycle_1
END SUBROUTINE SETNEB! ********** Print of Data on Meridian Plane
SUBROUTINE WRTF2(f,z1,isf)
DIMENSION f(0:20,0:6)
COMMON /IJKP/na,i1,i2,k1,jp,kp,jkp
CHARACTER*8 z1; CHARACTER*4 form(7),area(7)DATA form/’(1H ’,’2X I’,’2,4X’,’ 21(’,’’,’F7.3’,’)) ’/DATA area/’ -3P’,’ -2P’,’ -1P’,’’,’ 1P’,’ 2P’,’ 3P’/
OPEN(20,FILE=’OUTPUT.DAT’)３次元翼列ポテンシャル流れのプログラム
WRITE(20,60)z1,naform(5)=area(isf+4)
DO k=6,0,-1; WRITE(20,form)k,(f(i,k),i=0,20); ENDDO
WRITE(20,61)(i,i=0,20)
60 FORMAT(//1H 10X ’* * * ’, A8, ’ * * *’, 10X ’na =’, I3/)
61 FORMAT(1H 2X ’k/i=’, 21I7)
62 FORMAT(1H 2X I2, 4X 21F7.3)
END! ********** Print of 3D Data
SUBROUTINE WRTF3(f,z1,isf)
DIMENSION f(1029),ff(0:20,0:6)
COMMON /IJKP/na,i1,i2,k1,jp,kp,jkp
CHARACTER*8 z1; CHARACTER*4 form(7),area(7)DATA form/’(1H ’,’2X I’,’2,4X’,’ 21(’,’’,’F7.3’,’)) ’/DATA area/’ -3P’,’ -2P’,’ -1P’,’’,’ 1P’,’ 2P’,’ 3P’/
OPEN(20,FILE=’OUTPUT.DAT’)WRITE(20,60)z1,naform(5)=area(isf+4)
DO k=6,0,-1
FORALL(i=0:20,j=0:6)ff(i,j)=f(i*jkp+k*jp+j+1)
DO j=6,0,-1; WRITE(20,form)j,(ff(i,j),i=0,20); ENDDO
WRITE(20,61)(i,i=0,20)
WRITE(20,62)k
ENDDO
60 FORMAT(//1H 10X ’* * * ’, A8, ’ * * *’, 10X ’na =’, I3/)
61 FORMAT(1H 2X ’j/i=’, 21I7)
62 FORMAT(1H 60X ’k =’, I2/)
END
! ********** Solution of Simultaneous Linear Eqns Having Reduced Array by Gaussian Elimination
SUBROUTINE GAUSSB(a,b,n,m1,m2,det)
DIMENSION a(n,m1:m2),b(n)det=1.
cycle_1: DO k=1,n; i2=MIN0(n,k-m1); j2=MIN0(m2,n-k)
det=a(k,0)*detb(k)=b(k)/a(k,0); IF(k==n)CYCLE cycle_1
FORALL(j=1:j2)a(k,j)=a(k,j)/a(k,0)
DO i=k+1,i2; j0=k-ib(i)=b(i)-a(i,j0)*b(k)FORALL(j=1:j2)a(i,j+k-i)=a(i,j+k-i)-a(i,j0)*a(k,j)
ENDDO
ENDDO cycle_1
DO i=n-1,1,-1; j2=MIN0(m2,n-I)
DO j=1,j2b(i)=b(i)-a(i,j)*b(i+j)
ENDDO
ENDDO
END
! ********** Matrix Inversion by Gauss-Jordan Reduction
SUBROUTINE MATINV(a,aw,iw,n,det,e)
DIMENSION a(n,n),aw(n),iw(n)det=1.; e=AMAX1(e,1.E-5); FORALL(k=1:n)iw(k)=k
DO k=1,nw=a(k,k); IF(ABS(w)>e)
GOTO 21l=k; 22 l=l+1; IF(l>n)STOP 5555

57

58

非粘性流れの FEM 解析

w=a(l,k); IF(ABS(w)<=e)
GOTO 22
FORALL(j=1:n)
aw(j)=a(k,j); a(k,j)=a(l,j); a(l,j)=aw(j)ENDFORALL
iw0=iw(k); iw(k)=iw(l); iw(l)=iw0
det=-det; 21 det=w*det
a(k,k)=1.
FORALL(j=1:n)a(k,j)=a(k,j)/w
cycle_1: DO I=1,N; IF(i==k)CYCLE cycle_1
w=a(i,k); a(i,k)=0.
FORALL(j=1:n)a(i,j)=a(i,j)-w*a(k,j)
ENDDO cycle_1
ENDDO
cycle_2: DO j=1,n; l=j; IF(iw(j)==j)CYCLE cycle_2
23 l=l+1; IF(iw(l)/=j)
GOTO 23
FORALL(i=1:n)aw(i)=a(i,j); a(i,j)=a(i,l); a(i,l)=aw(i)ENDFORALL
iw(l)=iw(j)ENDDO cycle_2
END

このプログラムでは，翼列の下流領域のメッシュ分割と循環 ΓD の値，式 (11.61) の右辺 b の値が予測子
修正子法で決定されるため，全体の計算も反復する必要がある．この反復計算における解の収束状況を次
表に示す．
na

Drmax

Dthmax

DGmax

×104
DMmax

1
2
3
4
5
6
7
8
9
10

—
61
25
23
9
6
4
4
4
4

—
230
648
2892
964
354
580
101
239
94

—
4225
1677
2468
749
558
289
227
233
125

10152
1188
1679
449
182
161
69
68
42
38

Drmax, Dthmax, DGmax, DMmax はそれぞれ下流側の k 面の r，周期境界の θ，循環 ΓD ，流れ場のマッハ数 M
(n−1)
(n)
| で，表には 104 倍した値が示されている．このプログラムの計算は 10 回の最大変化量 max |fm −fmの反復でほぼ収束するが，それは解の収束状況を見ながらダンピング係数を調整したことによるもので，このような流れの計算では一般には更に多くの反復が必要である．ΓD の計算では強いダンピング係数 0.2 を課しているが，これは安定化のためというよりはその近くに解があるからである．これは，ΓD の値を少し変えると，下流側の分岐流線 (翼面上の流線に繋がる) も少し移動するが，このとき流線の方向が下流境界で少し変わるのに対し，後縁ではかなり大きく変わることによる．ダンピングをはずせば解は立ち所に発散する．Dthmax の値が na = 4 で急増するのはダンピング係数のかけ方による．
計算の結果は一応妥当なものと思われるが，上流と下流境界の近傍では多少ぎくしゃくする．それは羽根車から十分遠方とは言えない上流境界に一様な子午線速度と循環，下流境界に周方向に一様な流れを仮定
し，平衡条件の式 (11.63) が厳密に成立するとしていることによる．また φ の連立 1 次方程式の行列計算では多少の桁落ちが見られ，Dirichlet 条件が上流境界にのみ与えられるこの流れでは，下流境界の φ の値が−.006 程度小さくなる．ただしこの程度の桁落ちは，流速などの計算結果に全く影響しない．最大局所マッ

３次元翼列ポテンシャル流れのプログラム

59

ハ数は，負圧面側翼先端部の翼弦中央に現れ，その値は 0.988 になる．最後に計算に用いたメッシュと計算の結果を示す．上方のメッシュは１つの羽根車流路の正面図で 10 °上方側から見たものである．下方のメッシュは同じものの側面図で 20 °前方側から見たものである．前方のメッシュを濃く，翼表面のものを更に濃く，また上流と下流境界のものを赤みを帯びた色で示している．計算の結果は静圧分布，等マッハ数線，
相対流速ベクトルで示している．鳥瞰図 k = 0 (ハブ面), 3, 6 (ケーシング面) は翼間流路の計算結果，子午面図 j = 0 (圧力面＋周期境界), 3, 6 (負圧面＋周期境界) はハブ−ケーシング間流路の計算結果である．羽根車下流側の随伴渦面 (j = 0 と j = 6) を横切って静圧は等しく，流速ベクトルは大きさ等しく一般に方向の異なるものになる．なお子午面上のベクトルは相対流速の子午面成分であるが，回転流面上のベクトルは相対流速ベクトルを鳥瞰したもので，その大きさは相対流速の大きさに対応しない．

図 11.6: ３次元翼列流れの計算メッシュ，圧力分布と等マッハ数線