第13 章 非圧縮性流れの解法|MAC 型解法

３次元流れや乱流の計算には，流速と圧力の方程式を解く原始関数法 (primitive function method) が適している．この種の解法の基本は Navier-Stokes(NS) 方程式と連続方程式を解くものである．しかしながら，実際にこれらの方程式を解こうとすると，NS 方程式から流速を求めることはできるが，連続方程式から残りの変数，圧力を求めることができない．このことは，圧縮性流れの解法で連続方程式から密度，運動方程式から流速，エネルギー方程式から例えば岐点内部エネルギーをそれぞれ求めることができるのとは対照的である．このために非圧縮性流れの解法では何らかの戦略が必要になる．この章には最も一般的に用いられてきた MAC 型の解法について述べる．MAC 法では NS 方程式と，連続方程式に代わるものとして圧力方程式が解かれる．他に 非圧縮性流れの有力な解法としては，連続方程式に擬似圧縮項を加えた方程式

から圧力を計算する擬似圧縮性法がある．ところで NS 方程式と圧力方程式を解く方法では，解が連続の条件を満足しなくなる虞がある．MAC 法ではこの難点を，いわゆる食違い格子を用いまた圧力の差分方程式にある付加項を導入することによってみごとに克服している．なおこのアイディアは，曲線座標格子に拡張することはできるが，任意形状要素を用いる有限要素法に拡張することは難しい．そのために有限要素法では，NS 方程式と連続方程式をガレルキン法で同時に解く方法，ペナルティ関数法，ラグランジュ剰数法，ソレノイダルの重み関数を用いる方法など 各種の解法が考案されてきた．

13.1 MAC 法
Harlow-Welch によって提案された MAC(marker and cell) 法は1 ，はじめはラグランジュ座標系におけるマーカー粒子の動きを求めると同時にオイラー座標系における格子セルの方程式を解くものであった．その後マーカー粒子は自由表面にのみ置かれるようになり，通常は格子セルの方程式のみが解かれるようになった．ここには格子セルの方程式の解き方について述べる．MAC 法の基礎方程式は次の保存形で書かれた Navier-Stokes 方程式と圧力方程式である．

u
u
du
@u
+ r  uu = ;rp +  r2u +g
dt@t
r2 p = ;r  (r  uu )+ r  g(13.1)(13.2)

ただしこの式は無次元化されており，g は重力加速度などの流体に作用する物体力である．圧力方程式 (13.2)

は式 (13.1) の発散を取り，連続の条件

ru = 0(13.3)

を考慮することによって導かれたものである．MAC 法では定常流れも非定常流れも時間進行法で求められ，
式 (13.1) の初期値問題から u ，式 (13.2) の境界値問題から p が交互に繰り返し計算される．

1 Harlow, F.H. and Welch, J.E., Numerical calculaion of time dependent viscous incompressible ow of uid with freesurface, Phys. Fluids, Vol.8(1965), 2182{89.

12

非圧縮性流れの解法|MAC 型解法





; 01

4
y

6





U

u

-

; 00



; 10

4

X

x



P

p 

4

V

x y

v

図 13.1: MAC 法の格子計算には図 13.1 に示す食違い格子 (staggered grid) が用いられる．格子点上に x y ，格子セルの中心に

圧力 p，辺の中心に流速成分 u v が定義される．ここでは図に点線で囲んだ 4 点に同一添字を用いることとし，必要に応じ添字 X P U V を付けこれらの点を区別することにする．普通にセル中心に添字 ij ，辺の中心に i +1=2 j  i j +1=2  : : : を付ける場合，あるいはセル中心点に添字 P ，辺の中心に東西南北にちなんで n e w s，隣接セル中心に N E W S を付けることもある．式 (13.1) の差分方程式は，時間微分に Euler 前進差分，空間微分に 2 次中心差分を用いれば次のようになる．

hp00 ; p;10 in

x
h
i
n
n+1
n
v00
= v00
; t g00 + p00 ;yp0;1
un00+1 = un00 ; t f00 +(13.4a)(13.4b)

ただし

f00 =

(u2 )00 ; (u2);1 0 + (vu)01 ; (vu)00 ; u ; (g ) 
00
00
x
y
P

 P

X

X

x

(uv)10 ; (uv)00 + (v2 )00 ; (v2 )0 ;1 ; v ; (g )
00
00
x
y
u
; 2u00 + u10 + u0 ;1 ; 2u00 + u01
u = ;1 0

g00 =

X

P

X



P

y

00





x2

(13.4c)y 2

 は Laplace 差分演算子である．式 (13.4c) の u2  vu uv v 2 は 2 次中心差分では(u2 )00 = (u00 + u10)2=4  (v2 )00 = (v00 + v01)2=4 (vu)00 = (uv)00 = (u0 ;1 + u00)(v;1 0 + v00 )=4
P

X

P

X





のようになる．なおレ イノルズ数が小さくなく不安定になる場合には，2 次中心差分の代わりに後で示す上
流差分スキームを用いるべきである．

MAC 法の圧力の差分方程式は，直接 圧力方程式 (13.2) の差分近似を取るのではなく次のように導出される．
圧力方程式は Navier-Stokes 方程式の発散を取り連続の条件を考慮することによって導かれたが，圧力の差分方程式も，Navier-Stokes 差分方程式 (13.4) の発散を数値的に取って，すなわち f(13.4a)10;(13.4a)00 g=x +f(13.4b)01;(13.4b)00g=y と置いて導かれる．
n+1
n
D00
= D00
; t

h f10 ; f00

i; g00 ; t p+ g01y00

xu10 ; u00 v01 ; v00
D00 =
+ y
x

n

n

SMAC法3

n+1+1  v n+1 が連続の条件を満足するように 0 とする
また連続の条件については，D00
は次に求める流速 un
ij
ij

n
n
は本来は 0 になるべき量であるがすでに求めた流速 un vij
に対しては一般に 0 にはならないので
が，D00
ij

残しておく．結局圧力の差分方程式は次のようになる．

p00 = ;f10 ; f00 g01 ; g00 D00; y + tx(13.5)

なお MAC 法の差分方程式はもともとは有限体積法の定式化から導出されたものである．

MAC 法で解を時間ステップ t 進める計算手順は次のようになる．

 まず既知の u  p  f  g の値を用い式 (13.4) から陽的に u +1 の値を計算する．
 次にこの u +1 を用い f g D の値を更新する．
 それから式 (13.5) の連立 1 次方程式を反復法で解いて p の値を更新する．
n

n

n

n

n

n

次に MAC 法の特徴について述べる．食違い格子は，一見 解法を複雑にしているように思われるが，式(13.4) の圧力項と発散 D の表現を簡単にし ，結果として連続の条件 D00+1 = 0 を満足する圧力の差分方程n式をコンパクトにしている．この圧力の差分方程式は，普通の格子を用いた差分法や四辺形混合要素を用いた有限要素法に見られる spurious 誤差の発生を完全に防いでいる2 ．圧力の差分方程式 (13.5) と単に圧力方
程式 (13.2) を差分近似したものとを比較すれば，前者は後者に付加項 D00 =t+D00 を考慮したものであることが分かる．圧力の差分方程式は通常 反復法で解かれ，完全に収束しないうちに次のステップに移る．

このときの圧力の差分方程式 (13.5) の残差を r00 とすれば，次に求める u n+1 に対して連続の条件の誤差は
n+1
D00= ;t r00のように見積もられる．この誤差は良性のもので解が収束すれば 0 に近付く．一方 付加項のない圧力の差分方程式の場合には，u n+1 に対して連続の条件の誤差は
n+1
n
D00
= (1+ t )D00; t r00のようになる．この誤差は集積し 計算を不能にする．そのためにこの種の解法ではこの章の冒頭に述べた何らかの戦略が必須になる．
2 下図に食違い格子を用いる MAC 法，各要素の中心に p，4 頂点に u v の定義される四辺形混合要素を用いる有限要素法，通常格子を用いる差分法に対し ，連続の条件と圧力差分方程式を比較する．正方形要素の場合にはこの FEM の圧力差分方程式は+ ;1 1 + 1 ;1 + ;1 ;1 ; 4 00 = ; 10 ; 00 + 11 ; 01 ; 01 ; 00 + 11 ; 10 + 00
2 2
22のようになり，圧力 00 は隣の 10 ;1 0 01 0 ;1 とは無関係になり，斜め隣の 11と結び付く．この圧力の方程式は，+ が偶数の点のものと奇数の点のものが独立に解かれることになり，そのために圧力の値は checker-board 状に波打つ．長方形や任意四辺形要素でもこのような傾向は強く現れる．この波打ち誤差は spurious 誤差といわれ，波打ちを鎮めるために間に合せの手段として平滑化 (smoothing) が広く行われている．
p11

p

p

p

p

f

f

h

p

i

f

f

g

g

g

h

p

 p

 p

g

D

h

 p

t

p

 :::

pij

j

continuity
condition

u

r

r

r

r
v

v

rr

rr

rr

rr

r
u

r

r
r

u

v

b

pressuredierenceequationbpb@@ @;;@ ;;;; @;@; @@

b

b
b

b

b

b

b

p b

b

b

b

b

b
b

b

b

p

b
b
b

MAC methodFEMbordinary gridb4

非圧縮性流れの解法|MAC 型解法

13.2 SMAC 法
この節の説明の都合上，MAC 法は Navier-Stokes 方程式と圧力の Poisson 方程式

u n+1 = u n ; t (r uu + rp ;  r2u ;g )nr2 p = r u=t ;r ( r uu ;  r2u ;g )(13.6a)(13.6b)を用い，u n が与えられるときに式 (13.6b) から pn ，式 (13.6a) から u n+1 を次々に計算する陽解法 (explicit

method) と見ることにする．
Amsden-Harlow の SMAC 法 (simplied MAC method) は，Navier-Stokes 方程式 (13.6a) に時間分離法(time-splitting method) を適用し ，圧力そのものではなく圧力の増分 の Poisson 方程式を解くものである．その基礎方程式はu  = u n ; t (r uu + rp ;  r2u ;g )nu n+1 = u  ; t r

r2 = r u=t

(13.7a)(13.7b)(13.7c)pn+1 = pn +となる．計算には MAC 法と同様の食違い格子と差分近似式が用いられ，解を時間ステップ t 進める計算
手順は次の通りである．

 u が与えられるときに，式 (13.7a) を陽的に解いて中間の流速 u を計算
 式 (13.7c) を反復法で解いて圧力の増分 を計算，また p +1 を計算
 式 (13.7b) から u +1 を計算
nnn

式 (13.7a) と式 (13.7b) を加えれば

u n+1 = u n ; t f(r uu ;  r2u ;g )n + rpn+1 gとなるから，この SMAC 法では，Navier-Stokes 方程式の対流項と粘性項にオイラー前進差分，圧力項にオイ
ラー後退差分を適用していることになる．なお式 (13.7c) は式 (13.7b) の発散を取り，連続の条件r u +1 = 0nを課して導出されたものである．SMAC 法は MAC 法の特徴をそのまま受け継いでいる．すなわち食違い格子が用いられ，また圧力の差分方程式には一見 MAC 法のような付加項は見当たらないが，u を介して間接的に考慮されている．定常流れを時間進行法で解く場合に，解が収束すれば(r uu + rp ;  r2u ;g ) ! 0 nu ! u n r u  ! 0! 0

u n+1 ! u となり，式 (13.7a) { (13.7c) の差分近似式の残差はすべて 0 に近付き，修正値も 0 に近付くようになる．
ここで SMAC 法を用い，よく解かれているバックステップ流路 (backward-facing stepped duct) の流れを計算する．計算領域と格子を図 13.2 に示す．ステップの高さ 1 とし ，座標の原点をステップの角に取り，その格子番号を i j = 0 とする．計算格子は等間隔長方形食違い格子 (uniform rectangular staggered grid)

で，格子間隔は x = 1=3 y = 1=5 である．このとき上流と下流の境界は i0 =側下方境界と上方境界は j0 =;30 i1 = 90，また下流;5 j1 = 15 になる．上流境界に流速 u，下流境界に静圧 p，また流れのレ イ

ノルズ数 Re の値を与え，定常の層流として計算する．次に FORTRAN 95/90 free source form で書かれたプログラムを示す．

SMAC法

5

3
y

0
= ;1
x

= ;10

0

30

図 13.2: バックステップ流路の計算格子
PROGRAM MAIN!***********************************************************************************************! Flow Problem: Steady 2D incompressible Flow through Backward-Stepped Duct
! Numerical Scheme: SMAC, Rectangular Grid, Central-Differences or Chakravarthy-Osher TVD Scheme!***********************************************************************************************PARAMETER(i0=-30,i1=90,j0=-5,j1=15)
DIMENSION u(i0:i1,j0:j1),v(i0:i1,j0:j1),uX(i0:i1,j0:j1),vX(i0:i1,j0:j1),p(i0:i1,j0:j1), &phi(i0:i1,j0:j1),c(0:4,i0:i1,j0:j1),z(i0:i1,j0:j1),f(i0:i1,j0:j1),locf(2,8),fmax(8), &al(-2:8)
COMMON //na,Re,dx,dy,dt
CHARACTER*10 z1,z2(8)DATA naf,Re,dx,dy,dt/10000, 100., .33333, .2, .04/
!naf=na_max, Re=Reynolds number
OPEN(20,FILE='OUTPUT.dat')! initial values: parabolic velocity profile
DO j=0,j1-1 y=(j+.5)*dy!pustream side duct
FORALL(i=i0:0)u(i,j)=(2./3.)*y*(3.-y)
ENDDO
DO j=j0,j1-1 y=(j-j0+.5)*dy u1=(9./32.)*y*(4.-y)!downstream side duct
FORALL(i=1:i1)u(i,j)=((90.-i)*u(0,j)+i*u1)/90.
ENDDO! Recovery pressure of the diverging duct: dp=coef*(1-.75*.75)/2, far-upstream and -downstream! pressure gradient: -4/3/Re, -9/16/Re, pressure losses of the duct: Dp=(4/3*10+9/16*30)/Re-dp
CALL COEF(c,i0,i1,j0,j1)
na=0 100 na=na+1!na=iteration number
CALL COMPZ(u,v,z,i0,i1,j0,j1)!compute zeta
CALL EXPCD(u,v,uX,vX,p,z,i0,i1,j0,j1,locf,fmax)!compute u*,v* by explicit
CALL COMPP(u,v,p,phi,c,i0,i1,j0,j1,resP,locf,fmax)!compute phi,pCALL COMPU(u,v,phi,i0,i1,j0,j1,locf,fmax)!compute u,v! Record converging process and decide convergence
IF(na==1)WRITE(20,'(A46)')'
na
resNS
resP
outflow
CPU-time'IF(MOD(na,20)==0)THEN
i=i1 resNS=0.
resNS = AMAX1(fmax(1),-fmax(2),fmax(3),-fmax(4))
flow = 3.*dy/8.*(3.*u(i,-5)+u(i,-4)+(u(i,10)+3.*u(i,11) & !outlet flow rate
+3.*u(i,12)+u(i,13))+u(i,13)+3.*u(i,14))
DO j=-4,8,2 flow = flow+dy/3.*(u(i,j)+4.*u(i,j+1)+u(i,j+2)) ENDDOCALL CPU_TIME(sec)!CPU-time
WRITE(20,'(I5,4F10.5)')na,resNS,resP,flow,sec
umin=0. DO j=j0,j1 umin=AMIN1(umin,u(i1,j)) ENDDO!search outlet backward flow
IF(umin<0.)WRITE(20,'(A60)')'should extend the duct to remove outlet backward flow'
ENDIF
IF(resP>10.)
GOTO 110!divergence
IF(resNS<1.E-4.AND.resP<1.E-4)
GOTO 110!convergence
IF(na<naf)
GOTO 100!iteration! Output computational results
110 CONTINUE
DO k=1,46

非圧縮性流れの解法|MAC 型解法

SELECT CASE(k)CASE(1) z1='velocity u' FORALL(i=i0:i1,j=j0:j1)f(i,j)=uX(i,j)
CASE(2) z1='velocity v' FORALL(i=i0:i1,j=j0:j1)f(i,j)=vX(i,j)
CASE(3) z1=' pressure ' FORALL(i=i0:i1,j=j0:j1)f(i,j)= p(i,j)
CASE(4) z1='vorticity ' FORALL(i=i0:i1,j=j0:j1)f(i,j)= z(i,j)
ENDSELECTWRITE(20,'(/1H 10X A7,A10,A7,10X A4,I5/)')'***** ',z1,' *****','na =',na
DO ii=i0,i1-40,40
DO j=j1,j0,-1 WRITE(20,'(I3,41F8.4,I3)')j,(f(i,j),i=ii,ii+40),j ENDDO
WRITE(20,'(41I8/)')(i,i=ii,ii+40)
ENDDO
ENDDO
! Compute and output reattachment point(rap)irap=5 10 irap=irap+1 IF(u(irap,j0)<0.)
GOTO 10!irap: grid point near rap
ii=-3 12 ii=ii+1 i=irap+ii
du=u(i,j0) d2u=uX(i,j0+1)-2.*du!du=Delta u, d2u=Delta^2u
d3u=u(i,j0+1)-3.*uX(i,j0+1)+3.*du!d3u=Delta^3ua=2. 13 fa = du+(a-1.)/2.*(d2u+(a-2.)/3.*d3u)!fa=f(alpha) by cubic interpolation
IF(ABS(fa)<.0001)
GOTO 14
dfa = d2u/2.+(2.*a-3.)/6.*d3u!dfa=f'(alpha)a=a-fa/dfa
GOTO 13!by Newton method
14 al(ii)=a IF(a>0.14)
GOTO 12!al(iii): alpha of u=0
IF(al(ii)<0.01)ii=ii-1!irap+iii: adjacent grid point to rap
xrap = (irap+ii+al(ii)/(al(ii-1)-al(ii)))*dx!comp rap by linear extrapolation
WRITE(20,'(5X A7,F5.1)')'xrap = ',xrap! Output locations and values of maximum and minimum of residuals and divergence
DATA z2/'maxdu =','mindu =','maxdv =','mindv =','maxphi =','minphi =','maxdiv =','mindiv ='/FORALL(i=1:8)locf(1,i)=locf(1,i)-31 FORALL(i=1:8)locf(2,i)=locf(2,i)-6
DO i=1,8 WRITE(20,'(5X A10,2I4,3X E10.3)')z2(i),(locf(k,i),k=1,2),fmax(i) ENDDO
CLOSE(20)
CALL BSDuct_CG(uX,vX,p,z)
STOP
END PROGRAM MAIN
! ********** Compute u* and v* by explicit scheme using central-differences
SUBROUTINE EXPCD(u,v,uX,vX,p,z,i0,i1,j0,j1,locf,fmax)!PARAMETER(i0=-30,i1=90,j0=-5,j1=15)
DIMENSION u(i0:i1,j0:j1),v(i0:i1,j0:j1),uX(i0:i1,j0:j1),vX(i0:i1,j0:j1),p(i0:i1,j0:j1), &z(i0:i1,j0:j1),du(i0:i1,j0:j1),dv(i0:i1,j0:j1),w1(-30:90),lo(2),locf(2,8),fmax(8)COMMON //na,Re,dx,dy,dt
DO i=i0,i1 j00=0 IF(i>0)j00=j0
FORALL(j=j00+1:j1-1)uX(i,j)=(u(i,j-1)+u(i,j))/2.
!u_X(u at point X)
ENDDO
DO j=j0+1,j1-1 i00=0 IF(j>0)i00=i0
FORALL(i=i00+1:i1-1)vX(i,j)=(v(i-1,j)+v(i,j))/2.
!v_X
i=i1 vX(i,j)=v(i-1,j)!outletENDDOFORALL(i=i0:i1,j=j0:j1)du(i,j)=0. FORALL(i=i0:i1,j=j0:j1)dv(i,j)=0.
DO j=j0,j1-1 i00=0 IF(j>=0)i00=i0!computation of u*
FORALL(i=i00 :i1-1)w1(i)=(u(i,j)+u(i+1,j))*(u(i,j)+u(i+1,j))/4.
!w1=uu at P
FORALL(i=i00+1:i1-1)du(i,j)=(w1(i)-w1(i-1))/dx
i=i1 du(i,j)=u(i,j)*(u(i,j)-u(i-1,j))/dx!outletDO i=i00+1,i1
du(i,j)=du(i,j)+(vX(i,j+1)*uX(i,j+1)-vX(i,j)*uX(i,j))/dy &+(p(i,j)-p(i-1,j))/dx+(z(i,j+1)-z(i,j))/dy/Re!corrections of NS eq
ENDDO
ENDDO

SMAC法
DO j=j0+1,j1-1 i00=0 IF(j>0)i00=i0
FORALL(i=i00:i1-2)dv(i,j)=(uX(i+1,j)*vX(i+1,j)-uX(i,j)*vX(i,j))/dx
i=i1-1 dv(i,j)=(uX(i,j)+uX(i+1,j))/2.*(v(i,j)-v(i-1,j))/dx
ENDDO
DO i=i0,i1-1 j00=0 IF(i>=0)j00=j0
FORALL(j=j00:j1-1)w1(j)=(v(i,j)+v(i,j+1))*(v(i,j)+v(i,j+1))/4.
DO j=j00+1,j1-1
dv(i,j)=dv(i,j)+(w1(j)-w1(j-1))/dy&+(p(i,j)-p(i,j-1))/dy-(z(i+1,j)-z(i,j))/dx/Re
ENDDO
ENDDO
FORALL(i=i0+1:i1,j=j0 :j1-1)u(i,j)=u(i,j)-dt*du(i,j)
FORALL(i=i0:i1-1,j=j0+1:j1-1)v(i,j)=v(i,j)-dt*dv(i,j)
lo=MAXLOC(du) FORALL(k=1:2)locf(k,1)=lo(k) fmax(1)=MAXVAL(du)lo=MINLOC(du) FORALL(k=1:2)locf(k,2)=lo(k) fmax(2)=MINVAL(du)lo=MAXLOC(dv) FORALL(k=1:2)locf(k,3)=lo(k) fmax(3)=MAXVAL(dv)lo=MINLOC(dv) FORALL(k=1:2)locf(k,4)=lo(k) fmax(4)=MINVAL(dv)
END SUBROUTINE EXPCD! ********** Set up coefficients of Laplace difference operator
SUBROUTINE COEF(c,i0,i1,j0,j1)!PARAMETER(i0=-30,i1=90,j0=-5,j1=15)
DIMENSION c(0:4,i0:i1,j0:j1)
COMMON //na,Re,dx,dy,dt
rsdx=1./dx/dx rsdy=1./dy/dy
DO j=j0,j1-1 i00=0 IF(j>=0)i00=i0
FORALL(i=i00:i1-1)c(1,i,j)=rsdx c(2,i,j)=rsdx
c(3,i,j)=rsdy c(4,i,j)=rsdy c(0,i,j)=-2.*(rsdx+rsdy)
ENDFORALLi= i00 c(1,i,j)=0. c(0,i,j)=c(0,i,j)+rsdx
ENDDO
DO i=i0,i1-1 j=0 IF(i>=0)j=j0c(3,i,j)=0. c(0,i,j)=c(0,i,j)+rsdy
j=j1-1 c(4,i,j)=0. c(0,i,j)=c(0,i,j)+rsdy
ENDDO
END SUBROUTINE COEF

7
!computation of v*!near outlet!w4=vv at P

!corrections of NS eq!u* at U!v* at V!inlet & step, phi_x=0!bottoms, phi_y=0!top,
phi_y=0! ********** Solve bound value problem of Poisson eqn for pressure inclement phi
SUBROUTINE COMPP(u,v,p,phi,c,i0,i1,j0,j1,resp,locf,fmax)!PARAMETER(i0=-30,i1=90,j0=-5,j1=15)
DIMENSION u(i0:i1,j0:j1),v(i0:i1,j0:j1),p(i0:i1,j0:j1),phi(i0:i1,j0:j1),c(0:4,i0:i1,j0:j1), &a(120,3),b(120),rhs(-30:89,-5:14),w(j0:j1),lo(2),locf(2,8),fmax(8)COMMON //na,Re,dx,dy,dt
n1=120 n2=20 kf=4 al=1.0 bt=0.5!alpha:accel, beta:damping
FORALL(j=j0:j1-1)p(i1,j)=0.
! Dirichlet cond
! Compute rhs
DO i=i0,i1-1 j00=0 IF(i>=0)j00=j0 DO j=j00,j1-1phi(i,j)=0.
rhs(i,j)=((u(i+1,j)-u(i,j))/dx+(v(i,j+1)-v(i,j))/dy)/dt!rhs
ENDDO ENDDO! Solve simultaneous linear eqs by ADI-methodk=0 40 k=k+1!i-sweep
DO i=i0,i1-1 im1=MAX(i-1,i0) j00=0 IF(i>=0)j00=j0
DO j=j00,j1-1 k=j-j00+1!setup coefs & rhs of linear eqs
a(k,1)=c(3,i,j) a(k,2)=c(0,i,j)*al a(k,3)=c(4,i,j)b(k)=rhs(i,j)-c(1,i,j)*phi(im1,j)-c(2,i,j)*phi(i+1,j)-(1.-al)*c(0,i,j)*phi(i,j)

8

非圧縮性流れの解法|MAC 型解法

ENDDO
CALL GAUSS3(a,b,n1,j1-j00)!solve linear eqs
FORALL(j=j00:j1-1)phi(i,j)=b(j-j00+1)!phi at P
ENDDO
k=k+1
!j-sweep
DO j=j0,j1-1 jm1=MAX(j-1,j0) i00=0 IF(j>=0)i00=i0
DO i=i00,i1-1 k=i-i00+1!setup coefs & rhs of linear eqs
a(k,1)=c(1,i,j) a(k,2)=c(0,i,j)*al a(k,3)=c(2,i,j)b(k)=rhs(i,j)-c(3,i,j)*phi(i,jm1)-c(4,i,j)*phi(i,j+1)-(1.-al)*c(0,i,j)*phi(i,j)
ENDDO
CALL GAUSS3(a,b,n1,i1-i00)!solve linear eqs
FORALL(i=i00:i1-1)phi(i,j)=b(i-i00+1)!phi at P
ENDDO
IF(k<kf)
GOTO 40! Compute residuals of phi equations and determine pressure p
resp=0.
DO i=i0 ,i1-1 im1=MAX(i-1,i0) j00=0 IF(i>=0)j00=j0
DO j=j00,j1-1 jm1=MAX(j-1,j0)
res=c(1,i,j)*phi(im1,j)+c(2,i,j)*phi(i+1,j)+c(3,i,j)*phi(i,jm1) &+c(4,i,j)*phi(i,j+1)+c(0,i,j)*phi(i,j)-rhs(i,j)!residuals
resp=MAX(resp,ABS(res))!max residual of phi fdep(i,j)=p(i,j)+bt*phi(i,j)!corrected p
ENDDO ENDDO
lo=MAXLOC(phi) FORALL(k=1:2)locf(k,5)=lo(k) fmax(5)=MAXVAL(phi)lo=MINLOC(phi) FORALL(k=1:2)locf(k,6)=lo(k) fmax(6)=MINVAL(phi)END SUBROUTINE COMPP! ********** Compute u and v from u* and v*
SUBROUTINE COMPU(u,v,phi,i0,i1,j0,j1,locf,fmax)
DIMENSION u(i0:i1,j0:j1),v(i0:i1,j0:j1),phi(i0:i1,j0:j1),d(i0:i1,j0:j1),lo(2),locf(2,8),fmax(8)COMMON //na,Re,dx,dy,dt
DO j=j0,j1-1 i00=0 IF(j>=0)i00=i0
FORALL(i=i00+1:i1)u(i,j)=u(i,j)-dt*(phi(i,j)-phi(i-1,j))/dx!u at U
ENDDO
DO j=j0+1,j1-1 i00=0 IF(j>0)i00=i0
FORALL(i=i00:i1-1)v(i,j)=v(i,j)-dt*(phi(i,j)-phi(i,j-1))/dy!v at V
ENDDO
! Compute divergent=u_x+v_y at point P
FORALL(i=i0:i1,j=j0:j1)d(i,j)=0.
DO i=i0,i1-1 j00=0 IF(i>=0)j00=j0
FORALL(j=j00:j1-1)d(i,j)=(u(i+1,j)-u(i,j))/dx+(v(i,j+1)-v(i,j))/dy
ENDDO
lo=MAXLOC(d) FORALL(k=1:2)locf(k,7)=lo(k) fmax(7)=MAXVAL(d)lo=MINLOC(d) FORALL(k=1:2)locf(k,8)=lo(k) fmax(8)=MINVAL(d)END SUBROUTINE COMPU! ********** Compute vorticity
SUBROUTINE COMPZ(u,v,z,i0,i1,j0,j1)!PARAMETER(i0=-30,i1=90,j0=-5,j1=15)
DIMENSION u(i0:i1,j0:j1),v(i0:i1,j0:j1),z(i0:i1,j0:j1)
COMMON //na,Re,dx,dy,dt! Compute zeta=v_x-u_y at point X
FORALL(i=i0:i1,j=j0:j1)z(i,j)=0.
!zeta=v_x on inlet, outlet, top & bottoms
DO j=j0+1,j1-1 i00=0 IF(j>0)i00=i0
FORALL(i=i00+1:i1-1)z(i,j)=(v(i,j)-v(i-1,j))/dx!zeta=v_x
ENDDO
FORALL(j=j0+1:0)z(0,j)=(3.*v(0,j)-v(1,j)/3.)/dx!stepSMAC法9

DO j=j0+1,j1-1 i00=1 IF(j>0)i00=i0
FORALL(i=i00:i1)z(i,j)=z(i,j)-(u(i,j)-u(i,j-1))/dy
ENDDO
DO i=i0,i1 j=0 IF(i>0)j=j0z(i,j)=z(i,j)-(3.*u(i,j )-u(i,j+1)/3.)/dy
j=j1 z(i,j)=(3.*u(i,j-1)-u(i,j-2)/3.)/dy
ENDDO
END SUBROUTINE COMPZ!zeta=zeta-u_y!bottom walls containing corner point!top! ***** Solve simultaneous linear eqns with tri-diagonal matrix
SUBROUTINE GAUSS3(a,b,n1,n)DIMENSION 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,-1 b(k)=b(k)-a(k,3)*b(k+1) ENDDO
END

メインプログラムでは，はじめに計算に必要なデータ，初期・境界値が入力される．次に SMAC 法の反復計算がサブルーチンを引用するかたちで行われる．それから計算終了時点で流速 u v ，静圧 p などの計算結果が出力される．このプログラムでは，上流境界に平行平板間流路の十分発達した層流の流速分布，すなわ

ち放物線速度分布 (平均速度 1) が与えられる．このとき，後に述べる Courant 数 CFL = umax t=x < 1という条件は t < x=umax = 0:222 になる．しかしながら実際の計算では時間ステップ t は 0.04 程度以下でないと解は収束しない．また Reynolds 数 Re = 1= は，対流項を中心差分で近似するときには，文献等にも見られるように，約 70 以下でないと収束しない．ここでは t = 0:04 Re = 70 に選ぶことにする．
流れが平均流速でこの流路を通り抜ける時間は約 50 である．したがって最大反復数 naf は少なくも 1000のオーダー必要である．resNS = resNS は Navier-Stokes 差分方程式の最大残差値，resp = resP は圧力差分方程式の最大残差値を表すものである．flow = flow は Simpson の 1/3 と 3/8 公式を用いて計算した出口流量で本来は 3.0 になるべきものである．これらの量は解の収束状況を示す指標として，CPU 時間とともに出力される．収束判定は resNS < 0:0001 かつ resp < 0:0001 によって行われる．

再付着距離は u の出力をもとに作図で求めることもできるが，ここではメインプログラムの終わりに計算で求めることにする．再付着距離 xrap の位置は uij0 の値が負になる領域 5 < i < irap の下流側にある．まず配列 al を用意し ，各格子線 x = const: 上で u = 0 になる点を求める．下壁面からの無次元距離= y=(y=2) を導入すれば ，この点の位置 は次の 3 次方程式で表される．
h

i
1
1
u( ) = u0 + u0 + ( ; 1) 2 u0 + ( ; 2)3 u0 = 02
32
ただし u0 = 0 u0 = u 0   u0 = u 0 +1 ; 2u 0  3 u0 = u 0 +1 ; 3u 0 +1 +3u 0 ，添え字 i は省略．この式は 3 次のニュートン前進補間公式に基づくものである．この式で = 0 は下壁面上で流速 0 を意味し，
大括弧の中 = 0，すなわち

1
1f ( ) u0 + ( ; 1) 2 u0 + ( ; 2)3 u0 = 02
3は流れの中のある点で流速 0 を示す．この式の根 はニュートン法で容易に求めることができる．
+1 =; f ( )=f 0 ( )11
f 0 ( ) = 2 u0 + (2 ; 3)3 u0
2
6
j

n

n

n

j

n

X

j

j

j

X

j

10

非圧縮性流れの解法|MAC 型解法> 0 の範囲すなわち剥離泡内でこのようにしての値が求められる．再付着距離は= 0 になる x で次式から求められる．



xrap = irap + I +



dx

I

I

;1 ;

I

ただし irap + I は再付着点の上流側の隣接格子点の番号である．
以下には各サブルーチンを初心者向けにやや詳しく説明する．プログラム中のコメントも参照されたい．
プログラムを簡潔にするためにステップの上流部分と下流部分を分けずに同時に処理している．
サブルーチン COMPZ： 格子点 X 上の渦度
ijX

= (v ; u )
x

u ;1
v ;1
; u ;y= v ;xij

y ijXの値を計算．その計算式は領域の内点では

i

ij

j

(13.8)ij

また境界点 例えば下方壁面上の点では

0

i X; u 1 =3
= ;(u ) 0 = ; 3u 0 yi

(13.9)iy i X

この式は滑りなしの条件 u i0X = 0 を考慮して導かれたものである．なおこの渦度は Navier-Stokes 方程式rの粘性項  2 u =; ， r2 v = yxの計算に用いられる．

サブルーチン EXPCD： U 点の u と V 点の v  を式 (13.7a) から計算．配列 u(i,j) は u ，u または
u +1 のどれかである．また配列 du(i,j) = du は次式の大括弧の中の定常 Navier-Stokes 差分方程式の残差で，その値は解が収束した時点でほぼ 0 になる．式 (13.7a) の差分方程式は対流項を保存形の 2 次中心差
n
ij

n
ij

ij

ij

分で近似すれば次のようになる．

 1w ; w1uij = unij ; t

ij

i

x

;1 + w2 +1 ; w2 + p ; p ;1 + 
j

ij

ij

ij

i

j

ij

+1 ; ijyn(13.10a)

y
x
 3
n
3
4
4
w i+1j ; w ij w ij ; w ij ;1 pij ; pij ;1
i+1j ; ijnvij = vij ; t+
+ y ;  xx
y

(13.10b)

ただし w1ij = (u2 )ijP  w2ij = w3ij = (uv )ijX  w4ij = (v 2 )ijP である．また下流境界のところでは，反射しないように対流項を 1 次上流差分で近似すれば次のようになる．

uij = unij ; t uij(u ; u ;1 ) + w2 +1 ; w2 + p ; p ;1 + +1 ; ijyn(i = i1 )
x
y
x

n
(
vij ; vi;1j ) w4ij ; w4ij ;1 pij ; pij ;1
i+1j ; ijnvij = vij ; t uijV++ y ;  x(i = i1 ; 1)
x
y
ij

i

j

ij

ij

ij

i

j

ij

(13.11a)(13.11b)

上式では下流境界の静圧 pi1 j = 0 としている．MAXLOC(du) は duij の最大値の位置 (i j )，MAXVAL(du) は

duij の最大値を求める関数で，MINLOC(du) と MINVAL(du) も最小値に関する同様の関数である．
サブルーチン COMPP： の Poisson 方程式 (13.7c) の境界値問題を解き静圧の修正値 pn+1 を求める．この方程式は 2 階の楕円型で，したがって境界条件はすべての境界に関数またはその法線微分の値を与えれば必要十分である．式 (13.7c) の差分方程式は 2 次中心差分を用いれば次のようになる．

;1 ; 2 + +1 +

i

j

ij

x2

i

j

ij

またここでは下流境界に関数値
式 (13.7b) からx= 0 または;1 ; 2 +
ij

y 2

ij

 







+1 = 1 ui+1j ; uij + vij ;1 ; vij
t
x
y

(13.12)= 0，その他の境界に法線微分値を与えることにする．固定壁境界では= 0 となる．また上流境界でも u 0 = u0 = u 0+1 であるから = 0 と

i1 j
y

n
i j

i j

n
i j

x

陰解法化と Neumann 安定解析

11

なる． ij の連立 1 次方程式はここでは ADI 法で解かれる．この連立１次方程式の係数，Laplace 差分演算子の係数 c はサブルーチン COEF で設定され，ADI 法における３重対角行列の連立１次方程式はサブルーチ

ン GAUSS3 を引用し Gauss 消去法で解かれる．配列 a, b, c は３重対角行列の連立１次方程式の係数，右辺， ij の連立 1 次方程式の係数の記録に用いられる．AL =は ADI 法の加速パラメータ，BT =  は静圧

p の修正計算における減衰係数である．ADI 法による の計算は，解があまり収束しない段階で次の時間
ステップの計算に移るが，ここでは掃引数 kf = kf = 4 とし x 方向と y 方向を交互に各 2 回掃引している．
pn+1 = pn +  ．
サブルーチン COMPU： un+1 と v n+1 の修正値を式 (13.7b) から計算．式 (13.7b) の差分方程式は 2 次中心差分を用いれば次のようになる．

unij+1 = uij ; tij

; ;1 
i

j

x

n+1
vij
= vij ; t

ij

; ;1

(13.13)ijy

下表は解の収束状況を示したものである．表中の resN S は Navier-Stokes 差分方程式 (13.10) と (13.11) の;3 は流速 ui1 j から求めた出口流量
5の誤差である．この表中の数字はすべて 10 倍されている．解の収束を resN S < 0:0001 かつ resP < 0:0001
最大残差値, resP は圧力差分方程式 (13.12) の最大残差値, flow = flow

から判定することにすれば ，収束解は n = 3039 で得られることになる．流速は上流側流路中心の 1.5 が最大で，圧力は上流境界の 0.1541 が最大，ステップの角のところの;0:0838 が最小である．なお判定条件を

1 桁厳しくした場合には n = 4000 を少し超えるところで収束する．再付着距離 5.4 は変わらない．この計算は = 1:0， = 0:5 と置いて得られたものである． < 1:0 例えば 0.9 にとって ADI 法の計算を加速することは不安定になりできない．またこれらの係数の値を変えても収束にほとんど 影響しない．それはこの定常流れの収束が，圧力差分方程式の収束よりも，渦度が流れ場の隅々まで完全に拡散することによって完了するからであろう．
表 13.1: バックステップ流路流れの解の収束状況 (Re = 70，SMAC 法)
10;5nresN S

resP

100
1000
2000
3000
4000

2585
513
80
6
1

3176
298
15
13
1

flow

;1875
1020
346
285
309

平行平板間流路の圧力降下は上流側が 4=3Re = :01905，下流側が 9=16Re = :00803 である．計算の結果

はそれぞれ .0191, .0069 となる．下流側の数値が小さいのは，下流境界がステップの十分下流ではなく，まだ圧力回復効果が作用しているところにあるからである．すなわち出口の速度分布には流路拡大の効果が残っており，最大流速は 1.152 でその最終値 1.125 に向けてなお減速中である．

13.3

陰解法化と

Neumann 安定解析

MAC 法，SMAC 法などの陽解法 (explicit method) ではクーラン数 C = max(juj=x jvj=y)t の制約が厳しく，特に高レイノルズ数流れではこの制約は一層厳しいものになる．つまり時間間隔 t を大きく取ることができないので多くの反復計算を要することになる．これに対し陰解法 (implicit method) では，上流化，線形化，因子化を適切に行えば，クーラン数をかなり大きく取ることができる．SMAC 法でも，式 (13.7a)12

非圧縮性流れの解法|MAC 型解法

を陰解法のものに書き換えればクーラン数を大きく取ることができる．ここではまず Navier-Stokes 方程式と同様の性質を持つスカラー量の輸送方程式@ + f +d=0@t(13.14)r

の陰解法化について述べる．

式 (13.14) に梯形則 (trapezoidal law) を適用すれば次式が得られる．

n+1 ; nnn+1
(13.15)
t + (1 )( f + d) + ( f + d) = 0
この式は  = 0 ならば 陽解法，= 1=2 ならば Crank-Nicholson 法，= 1 ならば 完全陰解法になる．また
 = 1=2 のときにのみ時間に関して 2 次精度でほかの  では 1 次精度である3 ．式 (13.15) は，f = u と置くことができれば線形化 (linearization) することができ，更に n+1 = n +  n と置けば次の  形陰解法
;

r

r

(delta-form implicit method) の式に書き換えることがでる．
(1+ t  u) n = t( f + d)nr

;(13.16)rここで左辺の演算子 r は u だけでなく  n にも作用していることに注意すべきである．なおここでは d

は時間によらないものとし dn+1 = dn とした．定常流れを時間進行法で求める場合には，この式の右辺の括弧の中は解くべき方程式 (13.14) の残差，また  n は修正値を表す．反復計算で解が収束すればこの残差はゼロに近付き修正値もゼロに近付く．この式の右辺は解の精度に直接関係するので精度良く離散化しなければならないが，左辺の演算子は，解の精度に無関係で，安定化と計算量の軽減のために大胆に近似することが許される．通常左辺には，安定で簡単な 1 次上流差分が用いられ，拡散項 d は主要部のみが考慮され，また近似因子法が適用され１次元問題に分解して解かれる．

非定常流れの場合には，時間微分については，2 次精度の Crank-Nicholson 法や後退差分法 (Gear 法)，更には各種の Runge-Kutta 法，Adams-Bashuforth 法などが用いられる．また空間微分については，右辺だ
けでなく左辺の演算子も解の精度に直接影響するので，2 次精度以上のものを用いなければならない．それゆえ式 (13.16) を直接解くことはあまり得策ではない．ここではまず 2 次の Crank-Nicholson 法を基に，各時間ステップごとに Newton 法で反復計算する方法について述べる．式 (13.14) に Crank-Nicholson 法を適用すれば ，

n+1 ; nt

+ 12 ( f + d)n +( f + d)n+1 = 0
r



r

この式を Newton 法で反復計算するために

1+ 12 t u  (m) = ( (m;1)
r

;

;

n+1  (m) = (m;1) +  (m) と置けば次式が得られる．

n ) ; 1 t (r f + d)n +(r f + d)(m;1) 

2

ただし上添字 (m) は (n +1) 時間ステップの第 m 近似値を表し ， (0) =(13.17)n ，解が収束すれば n+1 = (m)となる．この式から得られる収束解の精度は左辺の演算子には無関係で，左辺は定常流れの場合と同様に大胆に近似することができる．
3 Taylor 展開 n+1 = n +t n +(1=2!)t2
t
を式 (13.15) に代入すれば，次式が得られる．

 f n+1 = f n +t(@ff=@t)n +(1=2!)t2 (@ 2f=@t2 )n +@@2
n+2(r f + d)]n + 3!1 t2 @t=0
2 t +3 (r f + d)] +@t t
この式に式 (13.14) の関係を用いれば， = 1=2 のときにのみ 2 次精度であることが分かる．

( t+

r

1
2!

ntt +(1=3!)t3 nttt +f + d)n + t

陰解法化と Neumann 安定解析

13

次に式 (13.14) に 2 次後退差分を用いれば ，n;1 ; 4 n +3 n+1

2t+ ( f + d)n+1 = 0r

上と同様に各時間ステップごとに Newton 法で反復計算することにすれば ，その計算式は次のようになる．


1+ 23 t u  (m) = 31 ( n;1 4 n +3 (m;1)) 23 t ( f + d)(m;1)



r

;

;

;

r

(13.18)この式から得られる収束解の精度も左辺の演算子には無関係で，左辺は定常流れの場合と同様に大胆に近似することができる．

ここでスカラー輸送方程式 (13.14) の解の性質について簡単に述べる．この方程式は線形化すれば

@@t + ur+d=0(13.19)となる．u は通常 流速で，演算子 @=@t+u r は時空間 x t における流跡線 dx=dt = u に沿う内微分 (interior

derivative) を表している．式 (13.14) は初期値問題として解かれるが，式 (13.19) を１つの流跡線に沿って
時間 t に関して積分すれば(x  t) = (x 0  t0 )

Zt

;

d dt

(13.20)t0となる．ただしこの式の積分路は流跡線に沿って取るものとする．d は一般に拡散項や発生項を意味するが，輸送方程式が同次で d = 0 の場合には，１つの流跡線に沿って(t) = (t0 ) となり，式 (13.19) の解は流跡線に沿って増減することなく速度 u で伝播する波になる．特に定常の場合には各流線に沿って定になる．今この波を局所的に平面波4で近似すれば次のように表される．

= f ( x vt)


;

の値は一

(13.21)

ただし f は波形， x ; vt は位相 (phase) と呼ばれ，また  は波面の単位法線ベクトル，v =  u は波面の法線速度である．式 (13.21) は d = 0 の同次方程式 (13.14) を満足する．なお一般に，微分方程式が同次の場合にはこのように基本解が任意関数で与えられる5 ．輸送方程式が拡散項を含み同次でない d 6= 0 の場合には，解(x  t) は，式 (13.20) から明らかなように，一つの流跡線に沿って速度 u で伝播する間に d の値により増減する波になる．しかしその基本的性質は同じということができる．

以上見てきた輸送方程式の初期値問題は，数学的には適正 (proper) で安定に解くことができる．しかしながらこれを離散近似した差分方程式の初期値問題は必ずしも安定に解くこはできず，その差分近似に際しては特別の注意が必要である．ここでは１次元の輸送方程式@ + @f =  @ + u @  = 0(13.22)@t @x @t @x
について考える．この方程式は，xt 空間内に 1 パラメータ族の流跡線 dx=dt = u を持ち6 ，各流跡線に沿って解の値が一定になることを示している．ここで一つの例を挙げれば ， は塩水濃度，式 (13.22) は管内の塩水の流れを表す式と見ることができる．管の断面積は変化してもよいが流れは各断面内で一様とする．
4 数学的表現である．一般に n 次元空間において面は n ; 1 次元，線は 1 次元である．なお 2 次元から n ; 2 次元の部分空間は多様体 (manifold) と呼ばれる．
5 各項の階数が同じ方程式を同次方程式 (homogeneous equation) という．例えば Laplace 方程式は 2 階の同次方程式である．
6 一般に F (x t ) = 0 は xt 空間内の１つの曲線を表すが，パラメータ  を変化させれば次々に曲線を発生させることができる．
このようにして作られる曲線の集合を 1 パラメータ曲線族 (one-parameter family of curves) という．F (x t   ) = 0 は 2 パラメータ曲線族の式である．

14

非圧縮性流れの解法|MAC 型解法

上流端から流入する塩水の濃度が変化すれば濃度分布は非定常になるが，この流れでは生成や拡散がない

ので，各流体粒子 (素片) の濃度

は 速度 u で輸送される間 一定に保たれる．

始めに離散近似した方程式の安定性を判別する von Neumann の方法 (von Neumann stability analysis)

について述べる．この安定解析では，解(x t) の Fourier 成分を

n = rn exp( i  ix)i(13.23)n = (xi  tn ) xi = ix tn = nt i = p;1 である．この波は波数  で，振幅iのように置く．ただし

関数 r の大きさが 1 より小ならば減衰し 1 より大ならば増幅する．差分方程式に式 (13.23) を代入すれば ，rとC

ut=x の関係式が得られる．C は安定性にかかわる重要なパラメータでクーラン数 (Courant
number) と呼ばれる．差分方程式が安定であるためには，すべての波数  に対し r < 1 であることが必要
で，そのための条件は r と C の関係式から導出される．
j j

ここで一例として，同次の式 (13.22) を Euler 前進法と 1 次上流差分で近似するときに得られる差分方

程式

n+1 ; n + C + ( n ; n ) + C ; ( n ; n ) = 0
i
i i;1
i+1 i
i

(13.24)の安定条件を導く．ただし C  = (C  jC j)=2．この式に式 (13.23) を代入すれば r と C の関係式

r 1+ C (1 cos x)+ C sin x = 0
;

j

j

;

i

が得られる．この式にその共役複素数を掛け実数化すれば，次の安定解の得られるための条件式が導かれる．

r2 = 1 + 2(C 2 C )(1 cos x) 1;j

これより，1 ; cos x

j

;



0 であるから，次の Neumann の安定条件が得られる．

C = u t=x 1

j

j

j

j

(13.25)

この条件は点 xi  tn+1 を通る流跡線が時間 t = tn に区間 xi;1  xi+1 ] を通らなければならないことを示している．時間間隔 t はこの条件を満足する範囲に取られるべきで，この条件が満足されないときには計算が不安定になる虞がある．
次に１次元スカラー輸送方程式に梯形則を適用し離散近似したもう少し一般的な差分方程式

n+1 + (f^n+1 ; f^n+1 ) = n ; (1 ; )(f^n ; f^n )iii+1=2 i;1=2
i+1=2 i;1=2を取り上げる．ただし(13.26)= x=t，また f^ は数値流束 (numerical ux) と呼ばれるもので，上の 1 次上流

差分の場合には対称性を持たせれば

f^i+1=2 = 21 (fi + fi+1 ) 12 ui+1=2  i+1=2;

j(13.27)ji+1 ; i ．これは良く知られた Roe スキーム7 の数値流束関数である．1 次上となる．ただし  i+1=2 =流差分の場合の式 (13.26) に式 (13.23) を代入すれば，r と C の関係式



r 1+  C sin x + C (1 cos x) = 1 (1 ) C sin x + C (1 cos x)
7

372.

fi

j

j

;

g

;

;

fi

j

j

;

g

Roe, P.L., Approximate Riemann solvers, parameter vectors, and dierence schemes, J. Comput. Phys., Vol.43, 357{

TVD 安定性と TVD スキーム15が得られる．この式は

C 1 (1 2) C (1 cos x)r2 = 1 21+2
 C (1+  C )(1 cos x)
j

;

jf

;
j

;

j

j

j

jg

j

;

;



1

のように書き換えることができ，右辺第 2 項の分母は常に正であるから，その分子も正でなければならな
いことになる．これより次の安定条件が導かれる．

(1 2) C
;

j

j 

1

(13.28)

この条件から，1 次上流差分の場合には， = 1 の完全陰解法は Courant 数 jC j によらず安定で，また  = 1=2の Crank-Nicholson 法は中立安定で Courant 数をかなり大きく取ることができることになる．
次に 2 次中心差分の場合には数値流束関数は

f^i+1=2 = ui+1=2 ( i + i+1 )=2となる．この場合には，r と C の関係式はr(1+ C sin x) = 1 (1 )C sin xi;i(13.29);となり，この式は

C 2 sin2 xr2 = 1 (2 1)
1+ 2C 2 sin2 x
;

;



1

のように書き換えられ，これより次の安定条件が導かれる．

 1=2(13.30)

この条件から 2 次中心差分の場合には，Courant 数 jC j によらず  = 0 のオイラー前進法では不安定，= 1=2の Crank-Nicholson 法では中立安定，= 1 の完全陰解法では安定である．しかしながら，前述の計算例が示

すように，Neumann 安定解析では不安定であるにもかかわらず，高くないレ イノルズ数 70，最大 Courant
数 0.24 では安定に解を求めることができた．また後述の例が示すように，Neumann 安定解析では無条件安定であっても実際の計算で不安定になることはよくあることである．

13.4 TVD 安定性と TVD スキーム
TVD 安定性と TVD スキームに関しては，第 8 章と第 16 章に詳しい説明があるので，ここにはその概要を述べるにとどめる．前節に述べたように，非定常１次元非圧縮性流れの輸送方程式 (13.22) は，xt 面上に
1 パラメータ族の流跡線を持ち，解 (x t) は各流跡線に沿って一定値を取る．総変化量 (total variation)TVはこの (x t) に対して次のように定義される．

TV (t)Zjx (x t)jdx

ただし積分区間は任意の１つの流跡線から他の流跡線まで取るものとする.
式 (13.31) の T V は例えば右図の (t) に対しては
次のようになる．

Z xb

j xjdxxa= (b ; 2 )+(1 ; min )+ jmin ; a j

TV =(13.31)16非圧縮性流れの解法|MAC 型解法

TV は一般に時間 t の関数であるがこの問題では時間によらず一定になる．TV はもともと非定常１次元圧
縮性流れのオイラー方程式に関して定義されたものである．この方程式の場合には，流跡線のほかに２つの1 パラメータ族の特性曲線が存在し，各特性曲線族に対する TV は時間によらず一定になる．しかし積分区
間に衝撃波がある場合には時間と共に減少する．すなわち TV は増加することはなく TVD(total variation
diminishing) である．TV の定義式 (13.31) は差分で表せば

TV ( n )X nn
j i+1 ; i ji(13.32)

また TVD 条件はTV ( n+1 ) TV ( n )(13.33)



となる．TVD スキームは常に TVD 条件 (13.33) を満足する解の得られる数値スキームのことである．

TVD 条件が満足されるときには，当然，(i) の新たな極値は発生しない，(ii) の最小値が減少，最大値が増加することもない．したがって TVD スキームを用いれば解を安定に求めることができる．
始めに 1 次上流差分の TVD 安定性について述べる．1 次元のスカラー輸送方程式 (13.26) に，1 次上流
差分の数値流束関数 (13.27) を代入し ，その両辺の TV を取れば ，3 条件

C + 0

C; 0 

(1 ) C



;

j

j 

1

(13.34)のもとで，次の関係が成立する．

TV ( n+1 ) TV (LHS ) = TV (RHS ) TV ( n )

これより 1 次上流差分は，式 (13.34) のもとで，TVD 条件 (13.33) を満足し ，TVD スキームであることが分かる．この場合の TVD 安定条件 (13.34) は，はじめの２つは常に満足され，(1;)jC j  1 にのみ留意す

れば良いことになる．例えば  = 1=2 の Crank-Nicholson 法の場合には，Neumann 安定条件ではクーラン
数 jC j を自由に選んでも良いということであったが，TVD 安定条件では jC j  2 となる．

2 次以上の差分スキームは，Courant 数 C を適切に選ぶだけでは不十分で不安定になることもある．以下には 2 次以上のスキームにおける解の不安定性と，それを TVD 安定化する対策について述べる．ChakravarthyOsher 型スキームは 2 次上流差分と 2 次中心差分を線形結合したもので，その数値流束関数は次のようにjj与えられる8 ．

1
1++
f^i+1=2 = f^i(1)+1=2 + 4 (1 )fi;1=2 + 4 (1+ )fi+1=2;

1 (1 )f ;1;(13.35)i+3=2 4 (1+ )fi+1=24
この式の右辺第 1 項は 1 次上流差分の数値流束で，第 2 項以降はこれを 2 次にする補正項である．また は結合のパラメータで， = 1 ならば 2 次上流差分，= 1 ならば 2 次中心差分，= 1=2 ならば Leonard のQUICK スキーム9 の数値流束になる．QUICK スキームの数値流束は上流側にシフトした 3 点を通る 2 次
多項式，すなわち ui+1=2 0 ならば点 fi;1  fi  fi+1 ，また ui+1=2 < 0 ならば点 fi  fi+1  fi+2 を通る 2 次
Lagrange 補間多項式の f (xi+1=2 ) である．式 (13.35) の数値流束を用いた差分式の精度は Taylor 展開をも
;

;

;

;

とに評価すれば次のようになる．

 @f 
1 (3 1)x2 @ 3f  + O(x3 )
^
f
)
=
+
x i+1=2 i;1=2@x i 12
@x3 i

1 (f^

;

;

(13.36)8 Chakravarthy, S.R. and Osher, S., A new class of high accuracy TVD schemes for hyperbolic conservation laws. AIAAPaper 85-0363, 1985.
9 Leonard, B.P., A stable and accurate convective modelling procedure based on quadratic upstream interpolation.
Computer Meth. Appl. Meth. Engng., Vol.19(1979),59{98.

TVD 安定性と TVD スキーム

17

最後の 2 項は打切り誤差で，これから  = 1=3 に取れば スキームの精度が 3 次になることが分かる．また

 が 1/3 に近い QUICK スキームは 2 次ではあるが，2 次中心差分や 2 次上流差分に比べ精度の高いものであることも分かる．

次にこれらのスキームの TVD 安定化について述べる．式 (13.35) の数値流束は 1 次上流差分の数値流束とこれを 2 次にする補正項からなる．1 次上流差分の部分は TVD 安定であるから，補正項の部分をこれらのスキームが TVD 不安定にならないようにその大きさを制限すれば，スキームは TVD になるはずである．
そのために各種の制限関数 (limiter) が用いられてきた．Chakravarthy-Osher TVD スキームでは minmod
制限関数が用いられその数値流束は次のように書かれる．

1
1
~+
+
f^i+1=2 = f^i(1)+1=2 + 4 (1 )f i;1=2 + 4 (1+ )f~i+1=2;;

1 (1 )f~;
i+3=2
4
;

;

1 (1+ )f~;i+1=2
4(13.37)ただし1f^i(1)+1=2 = 2 (fi + fi+1 fi+1=2 )
f~j+1=2 = minmod(fj+1=2  bfj;1=2 )
;j(13.38a)j(j = i i +1)(13.38b)f~j+1=2 = minmod(fj+1=2  bfj+3=2 ) (j = i 1 i)minmod(x y) = sign(x) max 0 min x  sign(x)y
また fj = ui+1=2 j  (j = i 1 i i +1 i +2) fi+1=2 = ui+1=2  i+1=2 である．
;

fj

;

j

(13.38c)

g

j

j

j

j

ここで式 (13.37) の数値流束の意味を考える．まず制限関数の働かない場合にはこの式はf^i+1=2 = fi+ + 41 (1 )(fi+ fi+;1) + 41 (1+ )(fi++1 fi+);;;+fi;+1 14 (1 )(fi;+2 fi;+1 ) 41 (1+ )(fi;+1 fi;)
(13.39)ただし fj = uj j  (j = i 1 i i +1 i +2) となる．すなわち一般には，f^i+1=2 の値は，u の正負により上下いずれかの式によって，上流側にシフトした 3 点 fi;1  fi  fi+1 または fi  fi+1  fi+2 の値から計算される．その値は u と  により，図 13.3 に示す範囲内にある．
しかしながら u の符号が u < 0 (x < xi + x) u 0 (xxi + x) のように変わる音速点では10 ，
;

;

;

;

;

;

=1
 = 1=3



fi;1

 fi
;
;

+2

;

;fi+1

HjJ^ @I  = ;1
fi

 fi
;
;

+2

=1



fi;1

;

;fi+1@RHY;; ;@I  = 1=3
fi ;  = ;1(u > 0)(u < 0)

図 13.3: f^i+1=2 の値の範囲
10 音速点で流速 u の符号が変わるのではなく，圧縮性 1 次元流れには３つの特性曲線族 dx=dt = u u  c が存在し ，その１つ圧力波の u ; c の符号が変わるということである．また一般に多次元流れでは，流速の横成分 v w の符号の変わるところで同様のこと
が起きる．

18

非圧縮性流れの解法|MAC 型解法

上式は次のようになる．

8>(;1 < < 0)
< fi + (1 ; )fi =4 + (1+ )(fi+1 ; fi )=4^
fi+1=2 = > (1+ )(fi+1 + fi)=4(0  < 1):
fi+1 + (1 ; )fi+1=4 ; (1+ )(fi+1 ; fi )=4 (1  < 2)

これらの式は  = 1 の場合には 2 次中心差分のままであるが， の値が 1 よりも小さくなると近似度は次第に低下する．例えば xi;1 < x < xi に音速点があり  = 1=3 の場合には，上式より 1 階微分の差分式は

1  1
1 1 
1 (f^
(fx)i = xi+1=2 f^i;1=2 ) = x 3 fi;1 + 2 fi + 3 fi+1;;となり，誤差評価するまでもなく破綻していることが分かる．念のため Taylor 展開をもとに誤差評価すれば ，この差分式は 2(fx )i =3+ fi =(2x)+    となり，誤差の大きさは O(x;1 ) である．2 次以上の高精度差分式を用いた計算結果が音速点の近傍で不自然に波打つことは良く知られている．
同様に u の符号が u0 (x < xi + x) u < 0 (x xi + x) のように変わるところでも
8>< fi+1 ; (1 ; )(fi+2 ; fi+1 + fi;1 )=4 ; (1+ )(fi+1 ; fi )=4 (;1 < < 0)^fi+1=2 = > (1 ; )(fi ; fi;1 ; fi+2 + fi+1 )=4 + (3 ; )(fi + fi+1 )=4(0  < 1):
fi + (1 ; )(fi ; fi;1 ; fi+2)=4 + (1+ )(fi+1 ; fi)=4(1  < 2)となり差分式は破綻する．これを回避するために各種の対策が立てられてきた．例えば 流束ベクトル分離法では，移行区間を設け，数値流束 f^; 100%の状態から，f^; の割合を徐々に減らし数値流束 f^+ の割合を

徐々に増やし ，f^+ 100%の状態に移行させるというものがる．この種の解法では不自然な波打ちは完全に回避できるが，その数値誤差は一般に小さくない．

Chakravarthy-Osher スキームの流束の式 (13.39) は，その右辺第 1 行または第 2 行それぞれが一種の補間式になっている．u の符号の変わるところでこれらの補間式が適当にミックスされたものは，上記のように補間機能を果たさない．補間機能を保つためには，２つの補間式のいずれか一方を選べば良く，ui+1=2 の符号によって選べば次式が得られる．

nf^i+1=2 = sign(u+i+1=2 ) fi + 41 (1 )(fi fi;1) + 14 (1+ )(fi+1 fi)
;

;

o

;

n
o
sign(u;i+1=2 ) fi+1 41 (1 )(fi+2 fi+1) 14 (1+ )(fi+1 fi)(13.40)u の符号の変わるところでこの式の精度を調べる．ui;1=2 < 0 ui+1=2 > 0 の場合には，この流束を用いた 1 階微分の差分式は (fi+1 fi;1 )=(2x) = (fx )i +O(x2 ) となる．また ui;1=2 > 0 ui+1=2 < 0 の場合には，1 階微分の差分式は (2 )(fi+1 fi;1 )=2x (1 )(fi+2 fi;2 )=4x = (fx )i +O(x2 ) となる．
いずれの場合にも  によらず 2 次精度が保証される．式 (13.40) は，式 (13.39) で fj = sign(u
i+1=2 ) fj 
(j = i 1 i i +1 i +2) と置いたものと同じである．式 (13.39) で fj = ui+1=2 j  (j = i 1 i i +1 i +2)と置けば ，もとの式と全く同じではないが，u の符号の変わるところでも  によらず 2 次精度の式が得られる．この Chakravarthy-Osher スキームの流束の式は  = 1=3 のときにも 2 次精度であるが，実際の計算の結果は 3 次精度の式を用いたものと違わない．それは u が準線形方程式 (13.22) の最高階数の微分ではな
;

;

;

;

;

;

;

;

;

;

;

;

j

;

j

;

くその係数であることを考えれば当然といえよう．

次に Chakravarthy-Osher TVD スキームの流束の式 (13.37) における制限関数の働きについて考える．

minmod 制限関数は式 (13.38c) で定義され，プログラムもこの式をもとに作られるが，分かり易く書き直せば次のようになる．

minmod(x y) =( x (x y 同符号 jxj  jy j)y (x y 同符号 x > y )0 (x y 異符号)
j

j

j j

TVD 安定性と TVD スキームb1

minmod(r b)minmod(r 1)

;

;
;
;

;

0

19

1 b

2

3

r

4

図 13.4: minmod 制限関数
また y=x = r と置けば

( 1 (rminmod(r 1) =1)

r (0 < r < 1)
0 (r 0)のようになる．この minmod(r 1) の値は図 13.4 に示すようになる．これより式 (13.38b) の意味するもの

は，当該勾配 fj +1=2 の大きさが隣の勾配の大きさを b 倍したものよりも小さいときにはそのまま，大きいときには bfj ;1=2 に置き換え，また隣の勾配の傾きが逆のときには 0 に置くというものである．このよ

うに勾配の大きさを制限するものを勾配制限関数 (slope limiter) という．b は 1 よりも大きい正数であるが，その値を大きく取れば制限関数が働き解の精度の低下する範囲は当然狭くなる．

次に 1 次元のスカラー輸送方程式 (13.26) に，Chakravarthy-Osher TVD スキームの数値流束関数 (13.37)を代入した式から求められる TVD 安定条件を結果のみ示す．

C~ ; 0 

C~ + 0(1 ) C~



;

j

j 

1

(13.41)

ただしC~  = u~t=x
C~ = C~ + C~ ;

jj(13.42a)(13.42b)

;

h
n u
1
;
i1=2 minmod(1 br ) ; minmod(r  b)o
u~i1=2  ui1=2 1 + 4
i1=2
i3=2
ui1=2oin u(13.42c)+ 1+4  i1=2 minmod(ri1=2  b) ; minmod(1 bri3=2 )ui1=2また rj+1=2 =  j +1=2 = j +1=21 である．これらの 3 条件の第 1 式と第 2 式より次の条件が導かれる．

1 b (3 )=(1 );(13.43);

また第 3 式より次の条件が導かれる．

(1 )b1 C
;

j

j 

1

b1 = 14 5  +(1+ )bf;g(13.44)

ただし jC j = (t=x) max(jui;1=2 j jui+1=2 j) である．Chakravarthy-Osher TVD スキームが文字通りのTVD スキームであるためには，式 (13.38b) の b が式 (13.43) を満足し ，また Courant 数 C すなわち時間
間隔 t が式 (13.44) を満足するように選ばれなければならない．

20

非圧縮性流れの解法|MAC 型解法

広く使われている 3 次上流差分 ( = 1=3) の場合には b  4 になる．その上限 b = 4 に取れば ，陽解法

( = 0) では C

0:4，Crank-Nicholson 法 ( = 1=2) では C 0:8 になり，完全陰解法 ( = 1) では Cの上限はない．また Crank-Nicholson 法で C = 1 に取れば，b 2:5 になる．3 次上流差分の場合に，b = 4に取れば ，式 (13.37) の数値流束は次のようになる．
8><fi + 1 minmod(fi;1=2  4fi+1=2 ) + 1 minmod(fi+1=2  4fi;1=2 ) (u  0)3^
fi+1=2 = > 6 1(13.45):fi+1 minmod(fi+3=2  4fi+1=2 ) 1 minmod(fi+1=2  4fi+3=2 ) (u < 0)63
次にこの式の制限関数の働きを図を使って説明する．ここでは u > 0，f > 0 の場合を取上げるがほかの場合も制限関数の働きは同じである．この場合には式 (13.45) の第 1 式は次のようになる．
nof^i+1=2 = fi + 61 minmod(r 4) + 13 minmod(1 4r) fi+1=2 (r = fi;1=2 =fi+1=2 )(13.46a)no= fi + 61 minmod(1 4r) + 13 minmod(r 4) fi;1=2 (r = fi+1=2 =fi;1=2 )(13.46b)
j

j 

j

j

;

j

j 

j

j



;

図 13.5 の左図は，fi  fi+1 を固定し fi;1 を変えた場合の f^i+1=2 の値を上式から計算したものである．r > 4

では第 1 の minmod 関数が働き，f^i+1=2 の値が fi+1 の値を超えないよう制限される．同様に右図は，fi;1  fi

を固定し fi+1 を変えた場合の f^i+1=2 の値を下式から計算したもので，r < 1=4 で第 1 の minmod 関数の働
き，f^i+1=2 の値が fi+1 の値を超えないように制限される．なおこれらの式の第 2 の minmod 関数はクーラン数の上限に関わるものである．

図 13.5: Chakravarthy-Osher TVD スキームの数値流束
ここで多少補足すれば，「 TVD スキームを用いた場合に最小値が減少，最大値が増加しない」というのは，1 次元の同次方程式 (13.22) の場合には輸送量 に対して，また１次元圧縮性流れの場合には３つの特性族の各特性に沿って伝播するエントロピーや２つのリーマン不変量に対して言っているのである．これらの不変量は，拡散項や発生項のある非同次の場合にも，拡散や発生の効果により増減し「 最小値が減少，最大値が増加しない」ということにはならない．また多次元の場合には，拡散や発生項がなくても，流れの x
成分には y z 成分が発生項の役割をするので同じことが言える．輸送方程式の実質微分の差分近似に TVD

スキームを用いれば ，拡散項や発生項の有無にかかわらず輸送量は良く保存され解は安定化する．TVD ス

キームを用いたために，本来最大値の増加する物理量が最大値が制限され増加しなくなるというようなことは有り得ないことである．

SMAC  形陰解法

2113.5 SMAC 形陰解法
SMAC  形陰解法は，SMAC 法の第 1 式 (13.7a) を式 (13.16) によって陰解法の式に書換えたもので，その全体の式は次のようになる．

1+ t ( un  2) uu = rhsn 
un+1 = u  t
2 n = u=t pn+1 = pn +
f

r

r

;

;

r

r

g

u = un + uur

(13.47a)(13.47b)(13.47c)u にも作用して
ただし rhs = ;t (r  uu + rp ;  r2u ;g ) である．式 (13.47a) の左辺の演算子 r は u

いることに注意すべきである． = 1=2 ならば時間に関し 2 次精度の Crank-Nicholson 法，= 1 ならば 1 次

精度のオイラー陰解法になる．なお式 (13.47b) と (13.47c) は SMAC 法のものと同じで，またこれらの式はSMAC 法と同じ手順で解かれる． 形陰解法の利点については 13.3 節を参照されたい．
次に非定常流れの SMAC  形陰解法の式は式 (13.17) を用い次のように書かれる．
1+ t ( un  2) uu(m) = (u(m;1) un )+ 12 (rhsn + rhs(m;1))
u (m) = u(m;1) + uu(m)(13.48a)u(m) = u(m) 21 t (m)(13.48b)2 (m) = 2 u(m)  p(m) = p(m;1) + (m)(13.48c)
t
f

r

;

;

r

r

g

;

;

r

r

この解法では非定常流れが 各時間ステップご とに反復計算によって求められ る．その予測値は u (0) =u n  rhs(0) = rhsn ，または時間に関する外挿によって求められる．

SMAC  形陰解法の式 (13.47a) または (13.48a) は近似因子法を適用することによって効率よく解くことができる．近似因子 (approximate-factorization, AF) 法は，旧ソ連の Yanenko によって考案された多次元問題を１次元問題に帰着させて解く部分段階法のひとつで，Holt(1971) 監修の訳本によって西側に紹介された11 ．一般に  形陰解法の式は，２次元の場合に，対流項に 1 次上流差分，拡散項に 2 次中心差分を用いれば形式的に次のように書くことができる．

(1+ c+1 x + c;1 x + c+2 y + c;2 y )u = rhsrr(13.49); + ;ただし c+
1  c1  c2  c2 は Courant 数に拡散係数を考慮したもの，rx  ry は後退差分，x  y は前進差分の+2演算子で，例えば c+
1 rx u = (C +=x )(u00 ;u;10) である．この式は差分式で表せば次のようになる．

c0 u00 c+1u;10 + c;1 u10 c+2u0;1 + c;2 u01 = rhs;;(13.50)

ただし c0 = (1+ jc1 j + jc2 j) である．式 (13.49) は従来の近似因子法では，(1+ c+1 x + c;1 x )(1+ c+2 y + c;2 y )u = rhsrr(13.51)のように因子化され，次式により２段階に分けて解かれる．

(1+ c1 )u00 c+1u;10 + c;1 u10 = rhs(1+ c2 )u00 c+2u0;1 + c;2 u01 = u00
j

j

;

j

j

;

(13.52a)(13.52b)

11 Yanenko, N.N., On the implicit dierence computing methods for solving multidimensional heat conduction equations (Russian), Izv. Uchebn. Zaved., Mathematica, 4(1963), 148{157. Yanenko, N.N., The Method of Fractional Steps.,
1971, Springer-Verlag.

22

非圧縮性流れの解法|MAC 型解法

この因子化に伴う誤差は，式 (13.52a) に式 (13.52b) を代入し因子化前の式 (13.50) と比較することによって次のように求められる．

c+2( c+1 u;1;1 + c1 u0;1 + c;1 u1;1) + c2 ( c+1 u;10 + c1 u00 + c;1 u10 )+c;2 ( c+1 u;11 + c1 u01 + c;1 u11 )
;

;

j

;

j

j

j

j ;

j

j

j

この誤差は c1  0 または c2  0 または u(x y ) が１次的に変化するところでは無視できるが ，一般には Courant 数の 2 乗の大きさになることが分かる．このため時間間隔 t を大きく取ると，解の収束はかえって悪くなり，さらに発散する虞もでてくる．この因子化の誤差を小さく抑えるために，Jameson らは

LU-SGS 法 (lower-upper symmetric-Gauss-Seidel method) 12 を提案している．この方法では式 (13.49) は(1+ c+1 x + c+2 y )c;0 1 (1+ c;1 x + c;2 y )u = rhsrrのように因子化され，次式により２段階に分けて解かれる．

c0 u00 c+1u;10 c+2u0;1 = rhs
c0 u00 + c;1 u10 + c;2 u01 = c0 u00;

(13.53a)(13.53b);

この因子化の誤差は，式 (13.53a) に式 (13.53b) を代入し因子化前の式 (13.50) と比較することによって次のように求められる．
;

c;0 1 c+1 c;2 u;11 +(c+1c;1 + c+2c;2 )u00 + c+2c;1 u1;1fg

この誤差はたかだか Courant 数の大きさであることが分かる．
式 (13.50) に対する近似因子法の式は c1  c2

d1 u00 d2 u;10 = d0 rhs ;

0 とすれば ，一般に次のように書くことができる．

d3 u00 d4 u0;1 = u00;

これら２つの式を１つに纏めた式と もとの式 (13.50) を比較すれば ，次の 4 つの条件が得られる．

d1 d3 = d0 c0  d2 d3 = d0 c1  d1 d4 = d0 c2  d2 d4 = 0
di はこれらすべての条件を満足させることはできないが，平均的に妥当なものとして次のように決めることができる．

d1 = d3 = d0 = c0  d2 = c1  d4 = c2
これより，c1  c2 < 0 の場合も含めて，筆者らの考案した次の改良型近似因子法 (IAF, improved approximate-factorization) の式を導出できる．

c0 u00 c+1u;10 + c;1 u10 = rhs
c0 u00 c+2u0;1 + c;2 u01 = c0 u00;;

(13.54a)(13.54b)

この因子化に伴う誤差は

c0 ;1 (c+1 c+2 u;1;1 c;1 c+2 u1;1 c+1 c;2 u;11 + c;1 c;2 u11 );;12 Jameson, A. and Yoon, S., Lower-upper implicit schemes with multiple grids for the Euler equations, AIAA J., 25(1987),
929{935. Yoon, S. and Kwak, D., Implicit Navier-Stokes solver for three-dimensional compressible ows, AIAA J., 30(1992),
2653{59.

SMAC  形陰解法

23

(1+ c+1 x + c;1 x + c+2 y + c;2 y + c+3 z + c;3 z )u = rhs(13.55)で Courant 数の大きさになる．

３次元の場合には  形陰解法の式はrrrとなる．従来の近似因子法の式は(1+ c1 )u000 c+1 u;100 + c;1 u100 = rhs(1+ c2 )u000 c+2 u0;10 + c;2 u010 = u000(1+ c3 )u000 c+3 u00;1 + c;3 u001 = u000
j

j

;

j

j

;

j

j

;

(13.56a)(13.56b)(13.56c)

LU-SGS の式は

c0 u000 c+1 u;100 c+2 u0;10 c+3u00;1 = rhs
c0 u000 + c;1 u100 + c;2 u010 + c;3 u001 = c0 u000;;(13.57a)(13.57b);

また改良型近似因子法の式はc0 u000 c+1 u;100 + c;1 u100 = rhs
c0 u000 c+2 u0;10 + c;2 u010 = c0 u000
c0 u000 c+3 u00;1 + c;3 u001 = c0 u000(13.58a)(13.58b)(13.58c);;;

ただし c0 = 1+ jc1 j + jc2 j + jc3 j である．近似因子法または改良型近似因子法の各式は 3 重対角行列の連立

1 次方程式で，Gauss 消去法で容易に解くことができる．ベクトルコンピュータでは，この計算は並べて係数を与え，並べて解くことによって並列処理すべきである．また LU-SGS 法では，２次元の場合には i + j
一定線上の点を並列処理すべきで，式 (13.53a) で i + j が増加する向きに，式 (13.53b) で i + j が減少する
向きに交互に掃引が行われる．また３次元の場合には i+j +k 一定面上の点を並列処理することになりかな

り煩わしい．これらの方法の因子化の誤差は次の表のようになる．
表 13.2: 因子化に伴う誤差の大きさ
２次元
近 似 因 子 法

LU-SGS 法

改良型近似因子法

c2
c
c

３次元

c3
c
c

ここで再び 13.2 節に示したバックステップ流路を通る定常２次元流れを SMAC  形陰解法で解くことにする．そのプログラムは，u  のサブルーチンを  形陰解法のものに換えたほかは前のものと同じである．
ただし陰解法では時間ステップを大きく取ることができるので t = 0:2 とし ，また安定性も良いのでサブルーチン COMPP でダンピングをはずし  = 1:0 とし ，漸近定常解が速く得られるようにしている．以下には  形陰解法のサブルーチン IMPCO のみ示すことにする．

! ********* Compute u* and v* by delta-form implicit method using Chakravarthy-Osher TVD scheme
SUBROUTINE IMPCO(u,v,uX,vX,p,z,i0,i1,j0,j1,locf,fmax)
DIMENSION u(i0:i1,j0:j1),v(i0:i1,j0:j1),uX(i0:i1,j0:j1),vX(i0:i1,j0:j1),p(i0:i1,j0:j1), &z(i0:i1,j0:j1),hf(i0:i1),ru(i0:i1,j0:j1),rv(i0:i1,j0:j1),rhs(2,i0:i1,j0:j1), &a(120,120,3),b(120,120),c(0:4,i0:i1,j0:j1),w1(i0:i1),w11(i0:i1),w2(j0:j1),&24

非圧縮性流れの解法|MAC 型解法

w21(j0:j1),lo(2),locf(2,8),fmax(8)COMMON //na,Re,dx,dy,dt
theta=1.
DO i=i0,i1 j00=0 IF(i>0)j00=j0
FORALL(j=j00:j1)w2(j)=u(i,j)
CALL INTUX(w2,w21,j0,j1,j00)
FORALL(j=j00:j1)uX(i,j)=w21(j)ENDDO
DO j=j0+1,j1-1 i00=0. IF(j>0)i00=i0
FORALL(i=i00:i1)w1(i)=v(i,j)
CALL INTVX(w1,w11,i0,i1,i00)
FORALL(i=i00:i1)vX(i,j)=w11(i)
ENDDO! Compute convection term using Chakravarthy-Osher TVD scheme
DO j=j0,j1-1 i00=0 IF(j>=0)i00=i0!hf=uu_P(uu at point P)
FORALL(i=i0:i1)w1(i)=u(i,j)
CALL INTUP(w1,w11,i0,i1,i00)!w11=u_P
i=i0
hf(i)=w11(i)*w11(i)!near inlet
i=0
uP=(12.*u(i+1,j)-u(i+2,j))/32. hf(i)=uP*uP!near step, u_0j=(u_x)_0j=0
i=i1-1 hf(i)=w11(i)*w11(i)!near outlet
i=i1
uP=(3.*u(i,j)-u(i-1,j))/2. hf(i)=uP*uP!nearby outside of outlet
cycle_1: DO i=i00+1,i1-1 uP=w11(i)m=1 IF(uP>0.)m=-1 ip=i+mIF(i==i1-1.AND.m==1)CYCLE cycle_1
du=u(i+1,j)-u(i,j) du1=u(ip+1,j)-u(ip,j)hf(i)=uP*uP+ABS(uP)*(-du/2.+AMINMOD(du1,4.*du)/6.+AMINMOD(du,4.*du1)/3.)
ENDDO cycle_1
FORALL(i=i00+1:i1)ru(i,j)=(hf(i)-hf(i-1))/dx!residuals of NSu eqn
ENDDO
DO i=i0+1,i1 j00=0 IF(i>0)j00=j0!hf=vu at X
j=j00
hf(j)=0.
!bottomj=j00+1 hf(j)=vX(i,j)*(3.*u(i,j-1)+2.*u(i,j)-u(i,j+1)/5.)/4. !bottom wall neighboring point
j=j1-1
hf(j)=vX(i,j)*(3.*u(i,j)+2.*u(i,j-1)-u(i,j-2)/5.)/4.
!neighbor top
j=j1
hf(j)=0.
!top
cycle_2: DO j=j00+1,j1-1
uX0=uX(i,j) vX0=vX(i,j)
m=1 IF(vX0>0.)m=-1 jp=j+mIF(J==j00+1.AND.m==-1 .OR. j==j1-1.AND.m==1)CYCLE cycle_2
du=u(i,j)-u(i,j-1) du1=u(i,jp)-u(i,jp-1)
hf(j)=vX0*uX0+ABS(vX0)*(-du/2.+AMINMOD(du1,4.*du)/6.+AMINMOD(du,4.*du1)/3.)
ENDDO cycle_2
FORALL(j=j00:j1-1)ru(i,j)=ru(i,j)+(hf(j+1)-hf(j))/dy!residuals of NSu eqn
ENDDO
DO j=j0+1,j1-1 i00=0 IF(j>0)i00=i0!hf=uv at X
i=i00
hf(i)=0.
!inlet and stepi=i0+1 hf(i)=uX(i,j)*(v(i-1,j)+v(i,j))/2.
!neighbor inlet
i=1
hf(i)=uX(i,j)*(3.*v(i-1,j)+2.*v(i,j)-v(i+1,j)/5.)/4.
!neighbor step
i=i1-1 hf(i)=uX(i,j)*(-v(i-2,j)+6.*v(i-1,j)+3.*v(i,j))/8.
!neighbor outlet
i=i1
hf(i)=uX(i,j)*(-v(i-2,j)+3.*v(i-1,j))/2.
!outletcycle_3: DO i=i00+2,i1-1 uX0=uX(i,j) vX0=vX(i,j)
m=1 IF(uX0>0.)m=-1 ip=i+mIF(i==i1-1.AND.m==1)CYCLE cycle_3
dv=v(i,j)-v(i-1,j) dv1=v(ip,j)-v(ip-1,j)hf(i)=uX0*vX0+ABS(uX0)*(-dv/2.+AMINMOD(dv1,4.*dv)/6.+AMINMOD(dv,4.*dv1)/3.)
ENDDO cycle_3
FORALL(i=i00:i1-1)rv(i,j)=(hf(i+1)-hf(i))/dx!residuals of NSv eqn
ENDDOSMAC  形陰解法
DO i=i0,i1-1 j00=0 IF(i>=0)j00=j0!hf=vv at P
FORALL(j=j0:j1)w2(j)=v(i,j)
CALL INTVP(w2,w21,j0,j1,j00)!w21=v_P
j=j00
hf(j)=w21(j)*w21(j)!near bottom
j=j1-1 hf(j)=w21(j)*w21(j)!near top
DO j=j00+1,j1-2 vP=w21(j)m=1 IF(vP>0.)m=-1 jp=j+m
dv=v(i,j+1)-v(i,j) dv1=v(i,jp+1)-v(i,jp)
hf(j)=vP*vP+ABS(vP)*(-dv/2.+AMINMOD(dv1,4.*dv)/6.+AMINMOD(dv,4.*dv1)/3.)
ENDDO
FORALL(j=j00+1:j1-1)rv(i,j)=rv(i,j)+(hf(j)-hf(j-1))/dy!residuals of NSv eqn
ENDDO! Compute pressure and viscous terms, and residuals of NS eqn
DO i=i0+1,i1 j00=0 IF(i>0)j00=j0
FORALL(j=j00:j1-1)ru(i,j)=ru(i,j)+(p(i,j)-p(i-1,j))/dx+(z(i,j+1)-z(i,j))/dy/Re
ENDDO
DO i=i0,i1-1 j00=0 IF(i>=0)j00=j0
FORALL(j=j00+1:j1-1)rv(i,j)=rv(i,j)+(p(i,j)-p(i,j-1))/dy-(z(i+1,j)-z(i,j))/dx/Re
ENDDO! Compute u* by implicit SMAC-scheme
ttx=theta*dt/dx tty=theta*dt/dy
FORALL(i=i0+1:0,j=j0:-1,l=1:4)c(l,i,j)=0.
FORALL(i=i0+1:0,j=j0:-1)c(0,i,j)=1.
DO i=i0+1,i1 j00=0 IF(i>0)j00=j0
FORALL(j=j00+1:j1-2)w2(j)=vX(i,j)
CALL INTVP(w2,w21,j0,j1,j00)!w21=v_U
DO j=j00,j1-1
up=(u(i,j)+ABS(u(i,j)))/2. um=(u(i,j)-ABS(u(i,j)))/2.
vp=(w21(j)+ABS(w21(j)))/2. vm=(w21(j)-ABS(w21(j)))/2.
c(1,i,j)=-ttx*(up+1./dx/Re)!coefs of linear eqs
c(2,i,j)= ttx*(um-1./dx/Re) IF(i==i1)c(2,i,j)=0.
!outletc(3,i,j)=-tty*(vp+1./dy/Re)c(4,i,j)= tty*(vm-1./dy/Re)c(0,i,j)=1.-c(1,i,j)-c(2,i,j)-c(3,i,j)-c(4,i,j)rhs(1,i,j)=-dt*ru(i,j)!rhs of linear eqs
ENDDO ENDDO! Computation du* by modified AF-method
DO j=j0,j1-1 l=j-j0+1 DO i=i0+1,i1 k=i-i0a(l,k,1)=c(1,i,j) a(l,k,2)=c(0,i,j) a(l,k,3)=c(2,i,j)b(l,k)=rhs(1,i,j)
ENDDO ENDDOCALL GAUSSZ(a,b,120,j1-j0,i1-i0)
FORALL(j=j0:j1-1,i=i0+1:i1)rhs(1,i,j)=b(j-j0+1,i-i0)!du**
DO j=j0,j1-1 l=j-j0+1 DO i=i0+1,i1 k=i-i0a(k,l,1)=c(3,i,j) a(k,l,2)=c(0,i,j) a(k,l,3)=c(4,i,j)b(k,l)=c(0,i,j)*rhs(1,i,j)
ENDDO ENDDOCALL GAUSSZ(a,b,120,i1-i0,j1-j0)
FORALL(j=j0:j1-1,i=i0+1:i1)rhs(1,i,j)=b(i-i0,j-j0+1)!du*
DO j=j0,j1-1 i00=0 IF(j>=0)i00=i0
FORALL(i=i00+1:i1)u(i,j)=u(i,j)+rhs(1,i,j)!u* at U
ENDDO! Compute v* by implicit SMAC-scheme
FORALL(i=i0:-1,j=j0+1:0,l=1:4)c(l,i,j)=0.
FORALL(i=i0:-1,j=j0+1:0)c(0,i,j)=1.
DO j=j0+1,j1-1 i00=0 IF(j>0)i00=i0
FORALL(i=i00:i1-1)w1(i)=uX(i,j)

25

26

非圧縮性流れの解法|MAC 型解法

CALL INTUP(w1,w11,i0,i1,i00)
i=0 IF(j<=0)w11(i)=(12.*uX(i+1,j)-uX(i+2,j))/32.
DO i=i00,i1-1
up=(w11(i)+ABS(w11(i)))/2. um=(w11(i)-ABS(w11(i)))/2.
vp=(v(i,j)+ABS(v(i,j)))/2. vm=(v(i,j)-ABS(v(i,j)))/2.
c(1,i,j)=-ttx*(up+1./dx/Re)
c(2,i,j)= ttx*(um-1./dx/Re) IF(i==i1-1)c(2,i,j)=0.
c(3,i,j)=-tty*(vp+1./dy/Re)c(4,i,j)= tty*(vm-1./dy/Re)c(0,i,j)=1.-c(1,i,j)-c(2,i,j)-c(3,i,j)-c(4,i,j)rhs(2,i,j)=-dt*rv(i,j)
ENDDO ENDDO! Computation dv* by modified AF-method
DO j=j0+1,j1-1 l=j-j0 DO i=i0,i1-1 k=i-i0+1a(l,k,1)=c(1,i,j) a(l,k,2)=c(0,i,j) a(l,k,3)=c(2,i,j)b(l,k)=rhs(2,i,j)
ENDDO ENDDOCALL GAUSSZ(a,b,120,j1-j0-1,i1-i0)
FORALL(j=j0+1:j1-1,i=i0:i1-1)rhs(2,i,j)=b(j-j0,i-i0+1)
DO j=j0+1,j1-1 l=j-j0 DO i=i0,i1-1 k=i-i0+1a(k,l,1)=c(3,i,j) a(k,l,2)=c(0,i,j) a(k,l,3)=c(4,i,j)b(k,l)=c(0,i,j)*rhs(2,i,j)
ENDDO ENDDOCALL GAUSSZ(a,b,120,i1-i0,j1-j0-1)
FORALL(j=j0+1:j1-1,i=i0:i1-1)rhs(2,i,j)=b(i-i0+1,j-j0)
DO j=j0+1,j1-1 i00=0 IF(j>0)i00=i0
FORALL(i=i00:i1-1)v(i,j)=v(i,j)+rhs(2,i,j)
ENDDO
lo=MAXLOC(uq) FORALL(k=1:2)locf(k,1)=lo(k) fmax(1)=MAXVAL(uq)lo=MINLOC(uq) FORALL(k=1:2)locf(k,2)=lo(k) fmax(2)=MINVAL(uq)lo=MAXLOC(vq) FORALL(k=1:2)locf(k,3)=lo(k) fmax(3)=MAXVAL(vq)lo=MINLOC(vq) FORALL(k=1:2)locf(k,4)=lo(k) fmax(4)=MINVAL(vq)
END SUBROUTINE IMPCO! ***** MINMOD function
FUNCTION AMINMOD(a,b)
s=SIGN(1.,a) AMINMOD=s*MAX(0.,MIN(ABS(a),s*b))
END FUNCTION! ***** Interpolate values at P from values at U
SUBROUTINE INTUP(u,u1,n0,n1,n00)
DIMENSION u(n0:n1),u1(n0:n1)
FORALL(i=n00+1:n1-2)u1(i)=(-u(i-1)+9.*(u(i)+u(i+1))-u(i+2))/16.
i=n00
u1(i)=(3.*u(i)+6.*u(i+1)-u(i+2))/8.
i=n1-1 u1(i)=(3.*u(i+1)+6.*u(i)-u(i-1))/8.
END SUBROUTINE INTUP! ***** Interpolate values at P from values at V
SUBROUTINE INTVP(u,u1,n0,n1,n00)
DIMENSION u(n0:n1),u1(n0:n1)
FORALL(j=n00+1:n1-2)u1(j)=(-u(j-1)+9.*(u(j)+u(j+1))-u(j+2))/16.
j=n00
u1(j)=(12.*u(j+1)-u(j+2))/32.
j=n1-1 u1(j)=(12.*u(j )-u(j-1))/32.
END SUBROUTINE INTVP! ***** Interpolate values at X from values at U
SUBROUTINE INTUX(u,u1,n0,n1,n00)
DIMENSION u(n0:n1),u1(n0:n1)
FORALL(j=n00+2:n1-2)u1(j)=(-u(j-2)+9.*(u(j-1)+u(j))-u(j+1))/16.
j=n00
u1(j)=0.

!w11=u_V!near step

!coefs of linear eqs!near outlet!rhs of linear eqs

!dv**

!dv*
!v* at V

SMAC  形陰解法

27

j=n00+1 u1(j)=(3.*u(j-1)+2.*u(j)-u(j+1)/5.)/4.
j=n1-1
u1(j)=(3.*u(j)+2.*u(j-1)-u(j-2)/5.)/4.
j=n1
u1(j)=0.
END SUBROUTINE INTUX! ***** Interpolate values at X from values at V
SUBROUTINE INTVX(u,u1,n0,n1,n00)
DIMENSION u(n0:n1),u1(n0:n1)
FORALL(i=n00+2:n1-2)u1(i)=(-u(i-2)+9.*(u(i-1)+u(i))-u(i+1))/16.
i=n00
u1(i)=0.
i=n00+1 u1(i)=(3.*u(i-1)+2.*u(i)-u(i+1)/5.)/4.
i=n1-1u1(i)=(3.*u(i)+6.*u(i-1)-u(i-2))/8.
i=n1
u1(i)=(3.*u(i-1)-u(i-2))/2.
END SUBROUTINE INTVX! ***** Solve simultaneous linear eqns with tri-diagonal matrix by Gauss elimination
SUBROUTINE GAUSSZ(a,b,m1,m,n)
DIMENSION a(m1,m1,3),b(m1,m1)
DO k=1,n-1 DO l=1,mb(l,k)=b(l,k)/a(l,k,2) a(l,k,3)=a(l,k,3)/a(l,k,2)b(l,k+1) =b(l,k+1) -a(l,k+1,1)*b(l,k)
a(l,k+1,2)=a(l,k+1,2)-a(l,k+1,1)*a(l,k,3)
ENDDO ENDDOFORALL(l=1:m)b(l,n)=b(l,n)/a(l,n,2)
DO k=n-1,1,-1
FORALL(l=1:m)b(l,k)=b(l,k)-a(l,k,3)*b(l,k+1)
ENDDO
END SUBROUTINE GAUSSZ

このサブルーチンでは，SMAC  形陰解法の式 (13.47a) の差分方程式がc1ij ui;1j + c2ij ui+1j + c3ij uij;1 + c4ij uij+1 + c0ij uij = rhs1ij
c1ij vi;1j + c2ij vi+1j + c3ij vij;1 + c4ij vij+1 + c0ij vij = rhs2ij
ucu i;1j + vcu ij ; vcu ij;1 + pij ; pi;1j +  ij+1 ; ij 
rhs1ij = ;t ucu ij ;xyxyj ; ucv ij + vcv ij ; vcv ij ;1 + pij ; pij ;1 ;  i+1j ; ij 
rhs2ij = ;t ucv i+1x
y
y
x

(13.59a)(13.59b)(13.59c)のように書かれ，左辺の演算子には 1 次上流差分，右辺の対流項には Chakravarthy-Osher TVD スキームが用いられる．すなわち  c2 = t u; ;  +
c1ij = t;u;ij x ij x
ij
x
x

  c4 = t v; ;   c0 = 1 ; c1 ; c2 ; c3 ; c4+;v;
c3ij = tij y ij yij
ij ij ij ij
ij y
y

(13.60a)n

ucu ij = (uijP )2 + juijP j ; 12 x uijP + 61 minmod(x ui;IjP  4x uijP )o+ 13 minmod(x uijP  4xui;IjP ) nvcu ij = vijX uijX + jvijX j ; 12 y uijX + 16 minmod(y uij;JX  4y uijX )o+ 13 minmod(y uijX  4y uij;JX )  : : : : : :(13.60b)28

非圧縮性流れの解法|MAC 型解法


ただし u
ij  vij は式 (13.59a) では U 点，(13.59b) では V 点の値が取られる．また ucu ij  : : : は Chakravarthy-

Osher TVD スキームの数値流束で，ucu ij  vcv ij は P 点上，vcu ij  ucv ij は X 点上に定義される．I = 1 (uijP 0) = ;1 (uijP > 0) J = 1 (vijX  0) = ;1 (vijX > 0) である．
このサブルーチンでは，右辺の対流項，圧力項，粘性項が計算され，連立 1 次方程式の係数と右辺の値が設定され，改良型近似因子法で u  v  の値が求められる．対流項の計算では境界条件を考慮する必要がある．ス;u2j )=32，
v1jX = (3v0j +2v1j ;v2j =5)=4 となる．下方と上方の境界に関しても同様である．改良型近似因子法は，ダクトの欠けている部分に 1 次式 1 u ij = 0 を補い，長方形領域に対して適用される．３重対角行列の連立 1

テップ (x = 0) のところの境界条件は u0j = v0j = (ux )0j = 0 で，3 次式を考えれば u0jP = (12u1j

次方程式はサブルーチン GAUSSZ を用いて解かれる．計算のベクトル化率を上げるために，次元の１つ多い

配列 a,b を用い，連立 1 次方程式の係数の設定と Gauss 消去法の計算はすべて並列に行われる．
解の収束状況を表 13.3 に示す．レ イノルズ数 Re = 70，時間間隔 t = 0:2，COMPP の= 1:0．また

resNS  resP  flow の定義は表 13.1 と同じで，これらの値はすべて 10は前と同じで，収束解は n = 706 で得られる．流速の最大値はダ クト入口中心の 1.5 で，圧力の最大値はダ クト入口の 0.1547，最小値はステップ角近傍の ;0:0840，再付着距離は 5.7 の近くにある．表 13.4 は時間間隔 t を変えた場合の解の収束に要する反復数 n を陽解法のものも含め比較したものである．この表か5 倍されている．解の収束判定条件

ら  形陰解法の反復数は陽解法に比べ 1/4 以下になることが分かる．なお CPU 時間は陽解法に較べそれほど 減らない． 形陰解法の利点は CPU 時間よりもむしろ高レ イノルズ数まで計算できることである．表

13.6 は高レ イノルズ数の計算結果をまとめたものである．更に高レ イノルズ数では剥離域が出口に達する

ため収束解は得られない．
表 13.3: バックステップダ クト流れの解の収束状況 (Re = 70，SMAC  形陰解法)
10;5

n
250
500
750
1000

resNS
508
62
6
1

resP
397
24
7
2

flow
66
243
362
293

表 13.4: バックステップダ クト流れの収束状況の比較 (Re = 70)
陽

解

法

形陰解法

Re
200
300
400
500

t
0.04
0.1
0.2
0.4
0.5
0.6
0.8

Cmax
0.18
0.45
0.9
1.8
2.25
2.7
3.6

n
3 034
1 171
706
655
700
1 106

pmax
0.15390.15410.15460.15440.15470.1540収束解得られず

pmin;0:0840;0:0849;0:0841;0:0849;0:0846;0:0849

表 13.5: バックステップダ クトの高レ イノルズ数流れ

n
1100
1884
2796
7854

再付着距離

11.7
15.3
18.4
20.8


1.0
1.0
1.5
1.5


1.0
0.5
0.5
0.5

備考

計算途中に出口逆流

直角座標系から一般曲線座標系への変換

29

13.6 直角座標系から一般曲線座標系への変換

直角座標 (cartesian coordinates) を x ，一般曲線座標 (general curvilinear coordinates) を で表す13 ．これらの座標間の変換の測度 (metric) x  : : : とx  : : : の間には次の関係が成立する14 ．
2
3 2
3 2
3
x x x
1 0 0
x y z
6y y y 7 67 6
74    5 4 x y z5 = 4 0 1 0 5z z z x y z
0 0 1(13.61)

また変換のヤコビアン J はx x xJ = y  y y z z z(13.62)となる．これらは２次元の場合には次のようになる．

x = J y  x = ;J y  y = ;J x  y = J x 

J = x y ; x y

x に関する微分を に関する微分に書き換えれば次のようになる15 ．
@ @j @@@xi = @xi @ j  r = (r j ) @ j(13.63)(13.64)

具体的には，例えばy@ =1 z@x J 
@=@

y
z
@=@

y
z 
@=@

２次元では

@ = 1 y @ ;y @ @x J  @  @

また一般に

@ @j
@ j J @xi = 0@
@ j Jr j = 0であるから次の関係が導かれる．

@ = J @ j @ = @ J @ j J @x
@xi @ j @ j @xii@
J r = @ J (r j )j
@
J r2 = @ Jgjk @@jk(13.65a)(13.65b)(13.65c)

ただし (g ij ) は測度テンソルで，その成分は次のように定義される．

@i@jgij = r i r j = @x@xkk

(13.66)13 ここでの一般は直交に対し 直交+非直交の意である．以前には非直交であることを強調し非直交曲線座標格子 (non-orthogonal
curvilinear coordinate grid) などと言われた．
14 例えば x = x(   ) = x( (x y z ) (x y z )  (x y z )) であるから，この式を y で微分すれば 0 = y x +y x +y x なる関係が得られ，式 (??) はこれらの関係を纏めたものである．
15 ここでは Einstein 総和規約 (Einstein cnovention) が用いられる．すなわち同じ 添字が繰り返されるときには ai bi = a1 b1 +a2 b2 + a3 b3 のように解釈する．

30

非圧縮性流れの解法|MAC 型解法

流線に沿う微分は次のように書き換えられる．

@ j @ = U @ = U  r~
u r = ui @x@ = ui @xj@@(13.67)

ただし Uj は反変速度 (contravariant velocity) のj 方向成分である．反変速度 U は 空間の流速で 16 ，物

i

i

j

j

理空間の流速 u との間には次の関係がある．

@ ju 
Uj = @xiiui = @x@ Uji(13.68)j

なおこのような関係を知らずに，基礎方程式の xi の微分を，式 (13.64) によって

j の微分に書き換えるとたいへん長い式になり収拾がつかなくなる．
流れの渦度  =ru の xi 方向成分はk
i = ijk @u@x(13.69)j

ただし (ijk ) は Eddington の  で8>< 1 ((ijk) = (123) (231) または (312))
ijk = > ;1 ((ijk) = (321) (213) または (132)):0(その他の場合)

例えば ，1 = @u3 =@x2

; @u2=@x3 である．任意のベクトル a に対し次の関係が成立する．

r l ra = r l r m  @@aa = J1 lmn @@xi @@ai = J1 lmn @ @ @@xx a
m
n m
m
n



(13.70)

なおこの式の導出には次式が用いられた．
2x = J y z      lmn @ @ @x = 0(13.71)z zm n
式 (13.70) で，特に a = u ならば ．次の反変渦度と反変速度の関係が得られる．
Zl = J1 lmn @ @ (gnj Uj )(13.72)m
ただし (gij ) は測度テンソルで，(gij );1 = (g ij )，その成分は次式で定義される．
k @xk
gij = @@xx  @@xx = @x
(13.73)@i j
i @ j
また Zl は反変渦度の l 方向成分で，反変渦度 Z と物理空間の渦度  の間には次の関係がある．
@ j    = @xi Z(13.74)
Zj = @x
i
i @ j
i
j

x +u

16 反変速度が写像空間の流速であることは次のように説明できる．ある時間に点 にある流体粒子は単位時間後に点に移動する．点 の曲線座標を とすれば点の曲線座標はr となる．したがって反変速度
r は写像空間で流体粒子が単位時間に移動する距離すなわち流速であることが分かる．

x



x +u

 +u 

x

U u 曲線座標格子の SMAC  形陰解法

3113.7 曲線座標格子の SMAC 形陰解法
物理空間の曲線座標格子を写像空間の立方体格子に写像する．立方体格子の稜の長さ      は通常 1 に取られるが，以下の説明では当面 省略せずに残しておくことにする．直角座標格子の MAC 型解法

を曲線座標格子の解法に拡張する際には，その特徴を損なわないようにしなければならない．そのために，各セル面の中心に妥当な流速成分として写像空間の流速である反変速度 U の成分またはこれに相当のものが定義される．既存の解法では，反変速度成分，体積流束，反変物理速度成分，共変物理速度成分が用いら~れているが，ここでは体積流束 (volume ux) U JUU を用いることにする．

連続方程式は式 (13.65) を用いれば次のようになる．

i = @ J @ j u  = @ U~ = 0
J @u@xi @ j @xi i @ j j

図 13.6 は(13.75)空間内の１つの立方体格子セルを示したものである．式 (13.75) は，この立方体格子セルにわたって積分すれば(U~100 ; U~000)  + (V~010 ; V~000 ) + (W~ 001 ; W~ 000 )  = 0(13.76)となる．式 (13.76) の中で，

u x x@xi 000U  =   v y yw z z 000U u000U  (x 010 ;x000 )  (x 001 ;x000 )

U~000   = J @ ui

;x 000 )(x 001 ;x 000 ) は x 空間内の対応する六面体格子セルの， = 0 面の面積ベクトルであ~るから，U000   はこの面を通る流量である．したがって式 (13.76) はこの六面体セルの６つの面から出入りする流量の和が 0 になるという意味の式である．写像空間内の流れに垂直に取られた単位面積当たり~ V~  W~ はこの体積流束の    方向成分で，の物理空間の体積流量は体積流束 (volume ux) と呼ばれ，U
また  =  =  = 1 ならば対応する六面体セルの面を通る流量になる．なお上式の導出には式 (??) か

ここに (x 010

ら得られる次の関係が用いられた．

1 y y    

x=J

z z

~l は流速 u の左から J
次に体積流束の運動方程式について述べる．体積流束成分 Ur l  を演算したもので~l の運動方程式は運動方程式 duu=dt  @uu=@t+r  uu = ;rp+ r2u +g の左から J r l  を演算あるから，U 6HHH
x001

HHH W~
HH


001

U~000HH 
HHj V~010:
6

:

:
U~100

V~000 HHHxj  H 6010 HHHjHx000 H
HHHW~ 000H
x100
図 13.6:

空間内の立方体格子セル

32

非圧縮性流れの解法|MAC 型解法

することによって導くことができる．その対流項は次のようになる17 ．

Cl = J r l  (ruu) = r l  @@ U~iu = @@ U~i Ul ; u  U~i @@ r l
i
i
i


@
= @ U~i Ul + ilj U~i Ujil

i j は Christoel の記号である．
@ 2x =  l  @xx   l  = 1 glk @gjk + @gik ; @gij 
ij 2@ i@ j i j @ l@i @j @kただし

(13.77)(13.78)なお式 (13.77) の導出には次の関係が用いられた18 ．

@xx  @ r = ;r  @ 2x = ;r   k  @xx = ; l 
l @ @l ij @ij@j @i li jk
粘性拡散項については 14.4 節に詳しく述べる．結局，体積流束成分 U~l の運動方程式は次のようになる．

@ U~ = ; @ U~ U ;  l U~ U ; g~li @p ;  @ (~g Z~ ) + J @ l g(13.79)
lij @
jk k
@t l
@ i i l ij i j
@i
@xk k
i
~i = JUi  g~ij = Jgij  g~ij = gij=J Z~ = JZ である．この式の右辺第 2 項は，対流項を保存形で表し
ただし Uたために出てきた付加項で，その値は滑らかな格子では通常小さく等間隔の直交格子では 0 になる．

13.2 節の SMAC 法の基礎方程式は一般曲線座標系の SMAC 法の式に書き換えれば次のようになる．

nU~l = U~ln ; t Cl +~gli @@p ; Dli@
U~ln+1 = U~l ; t g~li @i@ g~li @  = 1 @ U~   pn+1 = pn + @@t @ llil(l = 1 2 3)(13.80a)(13.80b)(13.80c)

SMAC 法の陰解法化は，前にも述べたように，NS 方程式の前段の式 (13.80a) だけを陰解法化すれば達成でき，それによって時間間隔 t を大きく取ることが可能になる．一般曲線座標系の SMAC  形陰解法の式は次のようになる．

noJ + t  U~i @@ ; Dl Ul = rhsnl (l = 1 2 3)(13.81a)irhsl = ;t Cl +~gli @@p ; Dl  U~l = U~ln + JUl
i@
U~ln+1 = U~l ; t g~li @(13.81b)i@ g~li @  = 1 @ U~   pn+1 = pn + (13.81c)@l @i
t @ l l~l = U~ln であるから，式 (13.81a) の右辺 rhs は残差，Ul は修正
定常流れの計算では，解が収束すれば U
値を表している．解が収束すればこの残差はゼロに近づき，修正値もゼロに近づく．したがってこの式の17 第 2 式から第 3 式へは式 (13.65) と (13.68)，第 3 式へは U~i Ul = rl (U~iu ) の微分が用いられた．

18 第 1 式から第 2 式へは @ @xx rl = @ jl = 0 なる関係，第 3 式へは上の Christoel 記号の式，それから第 4 式へは@i @jl (@xx=@k ) = kl を用いている．

r

@i

33

バックステップ流路流れのプログラム

右辺は，解の精度に直接関係するので精度良く離散化しなければならないが，左辺の演算子は大胆に近似することができ，1 次上流差分が用いられ，拡散項は主要部のみ考慮され，また近似因子法が適用される．

Dl は拡散項の近似演算子で無視することもできるが，岐点の近傍の流れではこの項の働きが大きいのでこ
6

こではその主要部のみ残すことにすれば ，一般に gii > gij (i = j ) であるから次のようになる．

n


o
D1 = e @@ g~22 g~33 @@ J  + @@ g~33 @@ g11 + @@ g~22 @@ g11(13.82)また２次元の場合にはno2D1 = e @@ g~22 @@  + @@ 2 (~g11)非定常流れの計算では， = 1=2 の Crank-Nicholson 法を用い，各時間ステップごとに Newton 反復法で反復計算を行うことにする．この非定常流れの SMAC  形陰解法の式は次のようになる．

n

o
;
 ;

J + t U~i @@ ; Dl Ul(m) = ; U~l(m;1) ; U~ln + 12 rhsnl + rhs(l m;1)  (l = 1 2 3) (13.83a)iU~l(m) = U~l(m;1) + JUl(m)(m)
U~l(m) = U~l(m) ; 21 t g~li @@i@ g~li @(m)  = 2 @ U~ (m) @@t @ llil(13.83b)p(m) = p(m;1) + (m)(13.83c)~l
上添字 (m) は時間ステップ (n +1) における第 m 近似値を意味する．U

(0)

(m)

~l
始め，解が収束すれば U= U~ln p(0) = pn と置いて計算を= U~ln+1  p(m) = pn+1 となる．一般に定常流れが収束するまでには数千回の反復が必要であるが，非定常流れの各時間ステップの解は 3  5 回の反復で十分である．定常流れの収束までに多くの反復を要するのは拡散効果によるもので，一方 1 時間ステップに拡散効果の実質及ぶ範囲は限られるので 1 時間ステップ当たりの反復回数は少しで済むことになる．

13.8 バックステップ流路流れのプログラム
前章と同じバックステップ流路を通る流れを，長方形格子の変わりに曲線座標格子を用い SMAC  形陰解法で解く１つのプログラムを示す．ただしレ イノルズ数がある程度高くなっても剥離域が下流境界に達しないように下流側流路はかなり延長されている．このプログラムは，長方形格子の SMAC  形陰解法のプログラムとだいたい同じであるがかなり長大になるので，始めに図 13.7 のフローチャートを使って計算手順の概要を説明する．このプログ ラムでは メインプ ログラムで Mode を指定することにより，定常ま

たは非定常流れ，低レ イノルズ数または高レ イノルズ数流れの都合４つの場合を選択できる．定常流れの初期値はこのプログラムの中で与えられるが，非定常流れの初期値は別途計算したものを読込ませなければならない．低レ イノルズ数流れは対流項を 2 次中心差分で計算するもの，高レ イノルズ数流れは 3 次のChakravarthy-Osher TVD スキームで計算するものである．

このフローチャートによれば，始めに GRID で格子形成，METRIC で後の計算に必要な測度，ヤコビアン，
測度テンソル，また COEF で圧力差分方程式の係数の設定が行われ，それから定常流れの予測値の計算または非定常流れの初期値の読込みが行われる．それから左側の反復計算に移るが，n は時間ステップ，m は非定常流れの各時間ステップにおける反復数である．各反復計算では式 (13.81) または (13.83) により COMPU1

34

非圧縮性流れの解法|MAC 型解法

~ ，COMPP で  と p，COMPU2 で U~ ，また COMPUZ で  の計算が行われる．n20 回ごとに NS 差分方程式
でUの最大残差値と圧力差分方程式の最大残差値を打出し ，また非定常流れでは CG に必要なデータをファイルに書込む．この反復計算は解が発散したときには直ちに中止し ，解が収束または規定の反復数に達した
時に次のステップまたは計算を終了する．終了直前には計算の最終結果をファイルに書込み CG で流れの可視化を行う．

START 
GRID

METRIC
COEF

PPP

P
Psteady?
PYes No
PREDCT  READ eCOMPUZn=0- n = n+1 m = 0-PPPPPPNo
undteady?

Yesm = m +1e

COMPU1COMPPCOMPU2COMPUZ

MOD(n,20)=0
 WRITE 
unsteady
 WRITE PPPYes
PPdiverge?
P
NoPPPPPPYesPP Yes
No 
PPsteady?
P
P


P

P
P No

? P

YesPPPm<mPNo eP?PP
Pn<n

P
Pe
Yes No WRITE
converge?
unsteady
f

f

CG


STOP 図 13.7: SMAC  形陰解法のフローチャート
以下にバックステップ流路を通る定常流れのプログラムを示す．
PROGRAM MAIN
! *****************************************************************************************! Flow Problem: 2D Incompressible Flow through Backward Facing Stepped Duct
! Numerical Method: General Curvilinear Coordinate Grid, SMAC Delta-Form Implicit Method,!
Chakravarthy-Osher TVD Scheme! *****************************************************************************************

バックステップ流路流れのプログラム

35

PARAMETER (if=140,jf=25)
DIMENSION x(2,0:if,0:jf),UJ(2,0:if,0:jf),u(2,0:if,0:jf),p(0:if,0:jf),phi(0:if,0:jf),&z(0:if,0:jf),UJX(2,0:if,0:jf),f(0:if,0:jf),co(0:if,0:jf,-1:1,-1:1),locf(2,8),fmax(8)COMMON //n,Re,dt,i0,j0,lowRe,steady,unsteady
COMMON /COMPUZ1/UJXLOGICAL lowRe,steady,unsteady
CHARACTER*4 z1,z2
CHARACTER*10 z3(8)DATA nf,mf,Re,dt,Mode/5000,5,1000.,1.,2/!nf=max n, mf=max(m), Re=Reynolds number
i0=if-1 j0=jf-1! Mode=1: for low Re steady flow
! Mode=2: for high Re steady flow
! Mode=3: for low Re unsteady flow
! Mode=4: for high Re unsteady flow
steady= (Mode==1 .OR. Mode==2)
unsteady = (Mode==3 .OR. Mode==4)
lowRe= (Mode==1 .OR. Mode==3)
CALL GRID(x,if,jf)!grid generation
CALL METRIC(x)!metrics, metric tensors
CALL COEF(co)!coefficients of press-fde
IF(steady)CALL PREDCT(x,UJ,u,p)!predict JU and p
IF(unsteady)READ(10)x,UJ,u,p!undteady: initial data
CALL COMPUZ(UJ,u,z,locf,fmax,CFLmax)!velocities vorticity & divergence
OPEN(20,FILE='OUTPUT.dat')
WRITE(20,60) IF(unsteady)WRITE(20,61)
n=0 100 n=n+1!n: step
m=0 101 IF(unsteady)m=m+1!m: approx
CALL COMPU1(m,UJ,u,p,z,locf,fmax)!compute JU^*CALL COMPP (UJ,p,phi,co,resP,locf,fmax)!compute phi and p
CALL COMPU2(UJ,phi)!compute JU
CALL COMPUZ(UJ,u,z,locf,fmax,CFLmax)!velocities vorticity & divergence! Decide convergence and output computational resultsIF(MOD(n,20)==0 .AND. (steady.OR.m==1))THEN
i=if resNS=0.
resNS=AMAX1(fmax(1),-fmax(2),fmax(3),-fmax(4))
flow=3./8.*(3.*UJ(1,i,0)+UJ(1,i,1)+UJ(1,i,23)+3.*UJ(1,i,24))!outlet flow rate
DO j=1,23,2 flow=flow+(UJ(1,i,j)+4.*UJ(1,i,j+1)+UJ(1,i,j+2))/3. ENDDOCALL CPU_TIME(sec)WRITE(20,62)n,resNS,resP,CFLmax,flow,sec!print residuals
ENDIF
60 FORMAT(/1H 3X 'n', 6X 'resNS', 6X 'resP', 7X 'CFLmax', 4X 'outflow', 5X 'CPU-time'/)61 FORMAT(1H 5X 'm', 6X 'resNS', 6X 'resP'/)
62 FORMAT(1H I5, 2F11.5, 2X F9.3, F11.5, 2X F9.2)
IF(resNS>2. .OR. resP>2.)
GOTO 110!diverge
IF(resNS<.0001 .AND. resP<.0001)THEN!converge
GOTO 110IF(steady)IF(unsteady)
GOTO 111
ENDIF
IF(unsteady.AND.m<mf)
GOTO 101!unsteady:(m)-approx
111 IF(n<nf)
GOTO 100!n-step
110 CONTINUE
DO k=1,6
SELECT CASE(k)CASE(1) FORALL(i=0:if,j=0:jf)f(i,j)=x(1,i,j)z1="
x" z2="" mz= 0CASE(2) FORALL(i=0:if,j=0:jf)f(i,j)=x(2,i,j)z1="
y" z2="" mz= 0CASE(3) FORALL(i=0:if,j=0:jf)f(i,j)=u(1,i,j)z1="
u" z2="" mz= 0CASE(4) FORALL(i=0:if,j=0:jf)f(i,j)=u(2,i,j)z1="
v" z2="" mz= 0CASE(5) FORALL(i=0:if,j=0:jf)f(i,j)=p(i,j)
z1="
p" z2="*10 " mz= 1

36

非圧縮性流れの解法|MAC 型解法

CASE(6) FORALL(i=0:if,j=0:jf)f(i,j)=z(i,j)
z1="zeta" z2="/10 " mz=-1
ENDSELECT!mz:scale factor
CALL WRTF(f,z1,z2,mz,n)!print out computational results
ENDDO
CALL REATTACH(x,UJ,UJX,if,jf,xrap,k)!compute reattachment point xrap
WRITE(20,'(5X A7,F5.1,3X I4)')'xrap = ',xrap,k!print out xrap
! Output locations and values of maximum and minimum of residuals and divergence
DATA z3/'maxUq =','minUq =','maxVq =','minVq =','maxphi =','minphi =','maxdiv =','mindiv ='/FORALL(k=1:2,i=1:8)locf(k,i)=locf(k,i)-1
DO i=1,8 WRITE(20,'(5X A10,2I4,3X E10.3)')z3(i),(locf(k,i),k=1,2),fmax(i) ENDDOCALL StepDF_CG(x,u,p,z,if,jf)!computer graphics for steady
CLOSE(20)
END PROGRAM MAIN
! ********** Generate curvilinear coordinate grid by analytical method
SUBROUTINE GRID(x,if,jf)!if=140,jf=25! Generate by sloving boundary value problem of elliptic pde
DIMENSION x(2,0:if,0:jf),x0(2,0:if,0:jf),xxi(4,0:if,0:jf),aJ(0:if,0:jf),g(3,0:if,0:jf), &P(0:if,0:jf),Q(0:if,0:jf),f(0:if,0:jf)
CHARACTER*4 z1
DATA dxi,det,al1/3*1./
nf=200 i1=35 i2=40 i0=if-1 j0=jf-1
OPEN(20,FILE='OUTPUT.dat')! ***** Give values of grid coordinates on boundaries by geometrical series! a+ar+ar^2+...+ar^(n-1)=a(r^n-1)/(r-1)
x(1,i1,0)=0. x(2,i1,0)=0.
DO i=i1-1,0,-1x(1,i,0)=x(1,i+1,0)-.09975*1.055**(i1-i-1) x(2,i,0)=0.
!bottom wall
ENDDO
DO i=i1+1,i2
x(1,i,0)=0. x(2,i,0)=x(2,i-1,0)-.10045*1.35**(i-i1-1)!step wall
ENDDO
DO i=i2+1,70x(1,i,0)=x(1,i-1,0)+.2*1.05**(i-i2-1)!bottom wall
ENDDO
FORALL(i=71:if)x(1,i,0)=x(1,70,0)+.8232*FLOAT(i-70)!bottom wall
FORALL(i=i2:if)x(2,i,0)=-1.
FORALL(i=0:if)x(2,i,jf)=3.
!top wall! ***** Give predict values of x on top bound
FORALL(i=0:i1)x(1,i,jf)=x(1,i,0)
FORALL(i=i1+1:i2)x(1,i,jf)=x(1,i1,jf)-x(2,i,0)
FORALL(i=i2+1:if)x(1,i,jf)=x(1,i2,jf)+(x(1,if,0)-1.)/x(1,if,0)*x(1,i,0)! ***** Determine starting values of x,y by interporation
DO i=0,if DO j=1,j0
et=j/FLOAT(jf) alp=(-et*et+1.6*et+.4)*et! The cubic polynomial is formed to satisfy the conds of f(0)=0, f(1)=1, f'(0)=.4, f'(1)=.6! for finer grid near walls. But it is impossible sufficiently to control grid spaces by it.
x(1,i,j)=(1.-alp)*x(1,i,0)+alp*x(1,i,jf)x(2,i,j)=(1.-alp)*x(2,i,0)+alp*x(2,i,jf)
ENDDO ENDDO
n=0 100 n=n+1 IF(n==10)al1=1.4! ***** Compute x_xi y_xi x_eta y_eta
DO l=1,2
FORALL(i=1:i0,j=0:jf)xxi(l ,i,j)=(x(l,i+1,j)-x(l,i-1,j))/2./dxi!x_xi, y_xi
FORALL(i=0:if,j=1:j0)xxi(l+2,i,j)=(x(l,i,j+1)-x(l,i,j-1))/2./det!x_eta, y_eta
FORALL(j=0:jf)xxi(l , 0,j)=-(3.*x(l, 0,j)-4.*x(l, 1,j)+x(l,2,j))/2./dxi
FORALL(j=0:jf)xxi(l ,if,j)= (3.*x(l,if,j)-4.*x(l,i0,j)+x(l,i0-1,j))/2./dxi

バックステップ流路流れのプログラム

37

FORALL(i=0:if)xxi(l+2,i, 0)=-(3.*x(l,i, 0)-4.*x(l,i, 1)+x(l,i,2))/2./det
FORALL(i=0:if)xxi(l+2,i,jf)= (3.*x(l,i,jf)-4.*x(l,i,j0)+x(l,i,j0-1))/2./det
ENDDO! ***** Compute J J^2 xi_x eta_x xi_y eta_y g_11 g_12 g_22
FORALL(i=0:if,j=0:jf)aJ(i,j)=xxi(1,i,j)*xxi(4,i,j)-xxi(3,i,j)*xxi(2,i,j)!Jg(1,i,j)=(xxi(1,i,j)*xxi(1,i,j)+xxi(2,i,j)*xxi(2,i,j))/aJ(i,j)!g_11/Jg(2,i,j)=(xxi(1,i,j)*xxi(3,i,j)+xxi(2,i,j)*xxi(4,i,j))/aJ(i,j)!g_12/Jg(3,i,j)=(xxi(3,i,j)*xxi(3,i,j)+xxi(4,i,j)*xxi(4,i,j))/aJ(i,j)!g_22/J
ENDFORALL! ***** Compute control functions P and Q! P: control i1-line and i2-line, Q: control grid spacing of j-direction
FORALL(i=0:if,j=0:jf)P(i,j)=0.
FORALL(i=0:if,j=0:jf)Q(i,j)=0.
FORALL(i=1:i1-1,j=1:jf)P(i,j)=2./(x(1,i+1,0)-x(1,i-1,0)) &!xi_xx*(1./(x(1,i+1,0)-x(1,i,0))-1./(x(1,i,0)-x(1,i-1,0)))
FORALL(i=i1+1:i2-1,j=1:15)P(i,j)=(15-j)/15.*2./(x(2,i+1,0)-x(2,i-1,0)) & !c*xi_yy*(1./(x(2,i+1,0)-x(2,i,0))-1./(x(2,i,0)-x(2,i-1,0)))
FORALL(i=i2+1:i0,j=1:10) P(i,j)=(10-j)/10.*2./(x(1,i+1,0)-x(1,i-1,0)) & !c*xi_xx*(1./(x(1,i+1,0)-x(1,i,0))-1./(x(1,i,0)-x(1,i-1,0)))
FORALL(j=1:jf)P(i1,j)=(P(i1-1,j)+P(i1+1,j))/2.
FORALL(j=1:jf)P(i2,j)=(P(i2-1,j)+P(i2+1,j))/2.
DO i=0,if DO j=1,jf
coef=25./(x(2,i,jf)-x(2,i,0))/(x(2,i,jf)-x(2,i,0))!coef=(d eta/d te)*(d ty/dy)^2
ty=(x(2,i,j)-x(2,i,0))/(x(2,i,jf)-x(2,i,0))!ty=(y-y_0)/(y_f-y_0)! Q(i,j)=coef*(6.*ty-4.)!te(ty)=ty^3-2ty^2+2ty
te'(0,.5,1)=(12/6 4.5/6 6/6)
Q(i,j)=coef*(8.*ty-5.)!te(ty)=(8ty^3-15ty^2+13ty)/6 ty'=(13/6
4/6 7/6)! Q(i,j)=coef*(10.*ty-6.)!te(ty)=(5ty^3-9ty^2+7ty)/3
te'=(14/6 3.5/6 8/6)
! Q(i,j)=coef*(12.*ty-7.)!te(ty)=(4ty^3-7ty^2+5ty)/2
te'=(15/6
3/6 9/6)! Q(i,j)=coef*(18.*ty-10.)!te(ty)=3ty^3-5ty^2+3ty
te'=(18/6 1.5/6 12/6)! for all te(ty), te(0, .5, 1.0)=(0, .625, 1.0),! Q=d^2eta/dy^2=(d ty/dy)^2*te''(ty)*(d eta/d te)=coef*te''(ty)! Lower Q's correspond to grids whose spaces are larger difference.
IF(i>i1.AND.i<i1+10.AND.j<10) &!modify grid spaces near corner
Q(i,j)=Q(i,j)-coef*12.*(1.-FLOAT(j)/10.)*(1.-ABS(i1+5-i)/5.)
ENDDO ENDDO! ***** Solve difference eqs of L(x)=-J(x_xiP+x_etaQ), L(y)=-J(y_xiP+y_etaQ) by SOR
FORALL(l=1:2,i=0:if,j=0:jf)x0(l,i,j)=x(l,i,j)
nn=0 110 nn=nn+1
DO 10 l=1,2 DO 10 i=1,i0 DO 10 j=1,j0
w=(g(3,i,j)*(x(l,i-1,j)+x(l,i+1,j))&-g(2,i,j)*(x(l,i-1,j-1)-x(l,i-1,j+1)-x(l,i+1,j-1)+x(l,i+1,j+1))/2. &+g(1,i,j)*(x(l,i,j-1)+x(l,i,j+1))&+aJ(i,j)*(xxi(l,i,j)*P(i,j)+xxi(l+2,i,j)*Q(i,j)))/2./(g(1,i,j)+g(3,i,j))
10 x(l,i,j)=x(l,i,j)+al1*(w-x(l,i,j))j=jf DO 11 i=1,i0!top wall
w=(g(3,i,j)*(x(1,i-1,j)+x(1,i+1,j))+g(1,i,j)*2.*x(1,i,j-1) &+aJ(i,j)*(xxi(1,i,j)*P(i,j)+xxi(3,i,j)*Q(i,j)))/2./(g(1,i,j)+g(3,i,j))
11 x(1,i,j)=x(1,i,j)+al1*(w-x(1,i,j))
i=0 DO 12 j=1,j0!inlet boundw=(g(3,i,j)*2.*x(2,1,j)+g(1,i,j)*(x(2,i,j-1)+x(2,i,j+1)) &+aJ(i,j)*(xxi(2,i,j)*P(i,j)+xxi(4,i,j)*Q(i,j)))/2./(g(1,i,j)+g(3,i,j))
12 x(2,i,j)=x(2,i,j)+al1*(w-x(2,i,j))!outlet boundi=if DO 13 j=1,j0
w=(g(3,i,j)*2.*x(2,i-1,j)+g(1,i,j)*(x(2,i,j-1)+x(2,i,j+1)) &+aJ(i,j)*(xxi(2,i,j)*P(i,j)+xxi(4,i,j)*Q(i,j)))/2./(g(1,i,j)+g(3,i,j))
13 x(2,i,j)=x(2,i,j)+al1*(w-x(2,i,j))

38

非圧縮性流れの解法|MAC 型解法

IF(nn<10)
GOTO 110! ***** Decide convergence of x,yadx=0.
DO i=0,if DO j=1,j0
adx=AMAX1(adx,ABS(x(1,i,j)-x0(1,i,j))+ABS(x(2,i,j)-x0(2,i,j)))ENDDO ENDDO
IF(n==1)WRITE(20,60)
IF(MOD(n,10)==0)WRITE(20,61) n,adx
60 FORMAT(/1H 3X 'n', 7X 'adx'/)
61 FORMAT(1H I5, F11.5)
IF(n<nf.AND.adx>1.E-4.AND.adx<5.)
GOTO 100ENDSUBROUTINE GRID!|dx|+|dy|! ********** Computation of metrics and metric tensors
SUBROUTINE METRIC(x)PARAMETER(if=140,jf=25,if1=280,jf1=50)
DIMENSION x(2,0:if,0:jf),x1(2,0:if1,0:jf1),xxi(4,0:if1,0:jf1),aJ(0:if1,0:jf1), &gU(3,0:if,0:jf),gV(3,0:if,0:jf),xix(4,0:if1,0:jf1),dxix(0:if,0:jf,2,4), &w1(0:if),w2(0:jf),w11(0:if1),w12(0:if1),w21(0:jf1),w22(0:jf1)
COMMON //na,Re,dt,i0,j0,lowRe,steady,unsteady
COMMON /METRIC1/x1 /METRIC2/aJ /METRIC3/gU,gV /METRIC4/xix /METRIC5/dxix! ***** Determine x1
DO l=1,2 DO j=0,jf
FORALL(i=0:if)w1(i)=x(l,i,j) CALL INTP(w1,w11,if)
FORALL(i=0:if1)x1(l,i,2*j)=w11(i)ENDDO ENDDO
x1(1,69,0)=(-x(1,33,0)+6.*x(1,34,0)+3.*x(1,35,0))/8. x1(2,69,0)=0.
x1(1,71,0)=0. x1(2,71,0)=(3.*x(2,35,0)+6.*x(2,36,0)-x(2,37,0))/8.
x1(1,79,0)=0. x1(2,79,0)=(-x(2,38,0)+6.*x(2,39,0)+3.*x(2,40,0))/8.
x1(1,81,0)=(3.*x(1,40,0)+6.*x(1,41,0)-x(1,42,0))/8. x1(2,81,0)=-1.
DO l=1,2 DO i=0,if1
FORALL(j=0:jf)w2(j)=x1(l,i,2*j) CALL INTP(w2,w21,jf)
FORALL(j=0:jf1)x1(l,i,j)=w21(j)ENDDO ENDDO! ***** Determine xxi, JDO l=1,2 DO j=0,jf1
FORALL(i=0:if1)w11(i)=x1(l,i,j) CALL DIFX(w11,w12,if1,.5)
FORALL(i=0:if1)xxi(l,i,j)=w12(i)!x_xi y_xi at point X
ENDDO ENDDO
DO l=1,2 DO i=0,if1
FORALL(j=0:jf1)w21(j)=x1(l,i,j) CALL DIFX(w21,w22,jf1,.5)
FORALL(j=0:jf1)xxi(l+2,i,j)=w22(j)!x_eta y_eta at X
ENDDO ENDDOFORALL(i=0:if1,j=0:jf1)aJ(i,j)=xxi(1,i,j)*xxi(4,i,j)-xxi(3,i,j)*xxi(2,i,j) !Jacobian J! ***** Determine g_11/J g_12/J g_22/J at U and V
DO i=0,if ii=2*i DO j=0,j0 jj=2*j+1
gU(1,i,j)=(xxi(1,ii,jj)*xxi(1,ii,jj)+xxi(2,ii,jj)*xxi(2,ii,jj))/aJ(ii,jj) !g_11/J at U
gU(2,i,j)=(xxi(1,ii,jj)*xxi(3,ii,jj)+xxi(2,ii,jj)*xxi(4,ii,jj))/aJ(ii,jj) !g_12/J at U
gU(3,i,j)=(xxi(3,ii,jj)*xxi(3,ii,jj)+xxi(4,ii,jj)*xxi(4,ii,jj))/aJ(ii,jj) !g_22/J at U
ENDDO ENDDO
DO i=0,i0 ii=2*i+1 DO j=0,jf jj=2*jgV(1,i,j)=(xxi(1,ii,jj)*xxi(1,ii,jj)+xxi(2,ii,jj)*xxi(2,ii,jj))/aJ(ii,jj) !g_11/J at VgV(2,i,j)=(xxi(1,ii,jj)*xxi(3,ii,jj)+xxi(2,ii,jj)*xxi(4,ii,jj))/aJ(ii,jj) !g_12/J at VgV(3,i,j)=(xxi(3,ii,jj)*xxi(3,ii,jj)+xxi(4,ii,jj)*xxi(4,ii,jj))/aJ(ii,jj) !g_22/J at V
ENDDO ENDDO! ***** Determine xix
DO i=0,if1 DO j=0,jf1

39

バックステップ流路流れのプログラム
xix(1,i,j)= xxi(4,i,j)/aJ(i,j) xix(2,i,j)=-xxi(3,i,j)/aJ(i,j)
xix(3,i,j)=-xxi(2,i,j)/aJ(i,j) xix(4,i,j)= xxi(1,i,j)/aJ(i,j)
ENDDO ENDDO! ***** Determine (xi_x)_xi
DO j=0,j0 jj=2*j+1
DO i=1,i0 ii=2*idxix(i,j,1,1)=xix(1,ii+1,jj)-xix(1,ii-1,jj)
dxix(i,j,1,3)=xix(2,ii+1,jj)-xix(2,ii-1,jj)
ENDDO
dxix(if,j,1,1)=3.*xix(1,if1,jj)-4.*xix(1,if1-1,jj)+xix(1,if1-2,jj)
dxix(if,j,1,3)=3.*xix(2,if1,jj)-4.*xix(2,if1-1,jj)+xix(2,if1-2,jj)
DO i=1,if ii=2*idxix(i,j,1,2)=xix(1,ii,jj+1)-xix(1,ii,jj-1)
dxix(i,j,1,4)=xix(2,ii,jj+1)-xix(2,ii,jj-1)
ENDDO ENDDO
DO i=0,i0 ii=2*i+1 DO j=1,j0 jj=2*jdxix(i,j,2,1)=xix(3,ii+1,jj)-xix(3,ii-1,jj)
dxix(i,j,2,3)=xix(4,ii+1,jj)-xix(4,ii-1,jj)
dxix(i,j,2,2)=xix(3,ii,jj+1)-xix(3,ii,jj-1)
dxix(i,j,2,4)=xix(4,ii,jj+1)-xix(4,ii,jj-1)
ENDDO ENDDO
ENDSUBROUTINE METRIC!xi_x, xi_y at X
!eta_x, eta_y at X!(xi_x)_xi!(xi_y)_xi

at U
at U

!(xi_x)_eta at U!(xi_y)_eta at U!(eta_x)_xi at V!(eta_y)_xi at V!(eat_x)_eta at V!(eta_y)_eta at V

! ********** Prediction of volume flux and pressure
SUBROUTINE PREDCT(x,UJ,u,p)
PARAMETER(if=140,jf=25,if1=280,jf1=50)
DIMENSION x(2,0:if,0:jf),UJ(2,0:if,0:jf),u(2,0:if,0:jf),p(0:if,0:jf),x1(2,0:if1,0:jf1), &aJ(0:if1,0:jf1),xix(4,0:if1,0:jf1)
COMMON //na,Re,dt,i0,j0,lowRe,steady,unsteady
COMMON /METRIC1/x1 /METRIC2/aJ /METRIC4/xix! ***** Set up starting values of u, U and p
DO j=0,jf ty=x1(2,0,2*j)/x1(2,0,jf1)
u(1,0,j)=6.*ty*(1.-ty)!mean velocity at inlet=1
ENDDO
DO j=0,j0 ty=x1(2,0,2*j+1)/x1(2,0,jf1)
UJ(1,0,j)=aJ(0,2*j+1)*xix(1,0,2*j+1)*6.*ty*(1.-ty)
FORALL(i=1:if)UJ(1,i,j)=UJ(1,0,j)
ENDDO
Dp=1./2.-.75*.75/2.
!recovery pressure in step ductp0=(4./3.*10.+9./16.*70.9)/Re-Dp!inlet pressure
FORALL(i=0:35,j=0:j0)p(i,j)=p0+4./3.*(x(1, 0,j)-x1(1,2*i+1,2*j+1))/Re
FORALL(i=100:i0,j=0:j0)p(i,j)=9./16.*(x(1,if,j)-x1(1,2*i+1,2*j+1))/Re
FORALL(i=36:99,j=0:j0)p(i,j)=(p(100,j)*(x1(1,2*i+1,2*j+1)-x1(1,71,2*j+1)) &+p(35,j)*(x1(1,201,2*j+1)-x1(1,2*i+1,2*j+1)))/(x1(1,201,2*j+1)-x1(1,71,2*j+1))
ENDSUBROUTINE PREDCT! ********** Compute JU^* and JV^* by delta-form implicit method using Chakravarthy-Osher!
type TVD scheme
SUBROUTINE COMPU1(ma,UJ,u,p,z,locf,fmax)
PARAMETER(if=140,jf=25,if1=280,jf1=50)
DIMENSION UJ(2,0:if,0:jf),Uq(2,0:if,0:jf),u(2,0:if,0:jf),p(0:if,0:jf),z(0:if,0:jf),aJ(0:if1,0:jf1),gU(3,0:if,0:jf),gV(3,0:if,0:jf),xix(4,0:if1,0:jf1),dxix(0:if,0:jf,2,4),
UJX(2,0:if,0:jf),hf(0:if),rhs(2,0:if,0:jf),rhs0(2,0:if,0:jf),UJ0(2,0:if,0:jf),a(140,140,3),b(140,140),c(0:4,0:if,0:jf),f(0:if,0:jf),w(0:if),w1(0:if),w11(0:if1),
w12(0:if1),w13(0:if1),w14(0:if),w2(0:jf),w21(0:jf1),w22(0:jf1),w23(0:jf1),w24(0:jf),
lo(2),locf(2,8),fmax(8)COMMON //na,Re,dt,i0,j0,lowRe,steady,unsteady

&
&
&
&
&

40

非圧縮性流れの解法|MAC 型解法

COMMON /METRIC2/aJ /METRIC3/gU,gV /METRIC4/xix /METRIC5/dxix /COMPUZ1/UJXLOGICAL lowRe,steady,unsteady
DATA theta,beta/1.,.5/!theta:trapezoidal, beta:damping
FORALL(l=1:2,i=0:if,j=0:jf)rhs(l,i,j)=0.
IF(steady.AND.lowRe)THEN!for training! ***** Compute convection term using central-differences
DO j=0,j0
FORALL(i=0:if)w1(i)=UJ(1,i,j) CALL INTP(w1,w11,if)!w11=JU_P(JU at poin P)
FORALL(i=0:i0)hf(i)=w11(2*i+1)*w11(2*i+1)/aJ(2*i+1,2*j+1)!hf=JUU_P
FORALL(i=1:i0)Uq(1,i,j)=hf(i)-hf(i-1)!JUU_xi
ENDDO
DO i=1,i0
FORALL(j=0:jf)hf(j)=UJX(2,i,j)*UJX(1,i,j)/aJ(2*i,2*j)!hf=JVU_X
FORALL(j=0:j0)Uq(1,i,j)=Uq(1,i,j)+hf(j+1)-hf(j)! +JVU_eta
ENDDO
DO j=1,j0
FORALL(i=1:i0)hf(i)=UJX(1,i,j)*UJX(2,i,j)/aJ(2*i,2*j)!hf=JUV_X
FORALL(i=1:i0-1)Uq(2,i,j)=hf(i+1)-hf(i)!JUV_xi
ENDDO
DO i=1,i0-1
FORALL(j=0:jf)w2(j)=UJ(2,i,j) CALL INTP(w2,w21,jf)!w21=JV_P
FORALL(j=0:j0)hf(j)=w21(2*j+1)*w21(2*j+1)/aJ(2*i+1,2*j+1)!hf=JVV_P
FORALL(j=1:j0)Uq(2,i,j)=Uq(2,i,j)+hf(j)-hf(j-1)! +JVV_eta
ENDDO
ELSE
! ***** Compute convection term using Chakravarthy-Osher TVD scheme
DO j=0,j0
FORALL(i=0:if)w(i)=UJ(1,i,j)/aJ(2*i,2*j+1)!w =U
FORALL(i=0:if)w1(i)=UJ(1,i,j) CALL INTP(w1,w11,if)!w11=JU_P
i=0
hf(i)=w11(2*i+1)*(w(i)+w(i+1))/2.
!near inlet
i=if-1 hf(i)=w11(2*i+1)*(w(i)+w(i+1))/2.
!near outlet
cycle_1: DO i=1,i0 UJP=w11(2*i+1)
m=1 IF(UJP>0.)m=-1 ip=i+mIF(i==i0.AND.m==1)CYCLE cycle_1
UP=(w(i+1)+w(i))/2. dU=w(i+1)-w(i) dU1=w(ip+1)-w(ip)hf(i)=UJP*UP+ABS(UJP)*(-dU/2.+AMINMOD(dU1,4.*dU)/6.+AMINMOD(dU,4.*dU1)/3.) !hf=JUU_P
ENDDO cycle_1
FORALL(i=1:i0)Uq(1,i,j)=hf(i)-hf(i-1)!JUU_xi
ENDDO! ***** Correct convection term just behind corner
FORALL(i=35:37)w(i)=UJ(1,i,0)*UJ(1,i,0)/aJ(2*i,1)
dum=w(36)-w(35) du=w(37)-w(36)
IF(UJ(1,36,0)<=0.)THEN Uq(1,36,0)=AMINMOD((dum/du+1.)/2.,2.)*du
ELSE Uq(1,36,0)=AMINMOD((du/dum+1.)/2.,2.)*dum ENDIF
DO i=1,i0
FORALL(j=0:j0)w(j)=UJ(1,i,j)/aJ(2*i,2*j+1)!w=U
j=0
hf(j)=0.
!bottom
j=1
hf(j)=UJX(2,i,j)*(3.*w(j-1)+2.*w(j)-w(j+1)/5.)/4.
!neighbor of bottom wall
j=j0 hf(j)=UJX(2,i,j)*(3.*w(j)+2.*w(j-1)-w(j-2)/5.)/4.
!neighbor top
j=jf hf(j)=0.
!top
cycle_2: DO j=1,j0 VJX=UJX(2,i,j)
m=1 IF(VJX>0.)m=-1 jp=j+m
IF(j==1.AND.m==-1 .OR. j==j0.AND.m==1)CYCLE cycle_2
UX=(w(j-1)+w(j))/2. dU=w(j)-w(j-1) dU1=w(jp)-w(jp-1)
hf(j)=VJX*UX+ABS(VJX)*(-dU/2.+AMINMOD(dU1,4.*dU)/6.+AMINMOD(dU,4.*dU1)/3.) !hf=JVU_X
ENDDO cycle_2
FORALL(j=0:j0)Uq(1,i,j)=Uq(1,i,j)+hf(j+1)-hf(j)! +JVU_etaバックステップ流路流れのプログラム
ENDDO
DO j=1,j0
FORALL(i=0:i0)w(i)=UJ(2,i,j)/aJ(2*i+1,2*j)!w=V
i=0hf(i)=0.
i=1hf(i)=UJX(1,1,j)*(w(i-1)+w(i))/2.
!neighbor inlet
i=i0 hf(i)=UJX(1,i,j)*(-w(i-2)+6.*w(i-1)+3.*w(i))/8.
!neighbor outlet
cycle_3: DO i=2,i0 UJX0=UJX(1,i,j)
m=1 IF(UJX0>0.)m=-1 ip=i+mIF(i==i0.AND.m==1)CYCLE cycle_3
VX=(w(i-1)+w(i))/2. dV=w(i)-w(i-1) dV1=w(ip)-w(ip-1)
hf(i)=UJX0*VX+ABS(UJX0)*(-dV/2.+AMINMOD(dV1,4.*dV)/6.+AMINMOD(dV,4.*dV1)/3.) !hf=JUV_X
ENDDO cycle_3
FORALL(i=0:i0-1)Uq(2,i,j)=hf(i+1)-hf(i)!JUV_xi
ENDDO
DO i=0,i0-1
FORALL(j=0:jf)w(j)=UJ(2,i,j)/aJ(2*i+1,2*j)!w =V
FORALL(j=0:jf)w2(j)=UJ(2,i,j) CALL INTVP(w2,w21,jf)!w21=JV_P
j=0
hf(j)=w21(2*j+1)*(12.*w(1)-w(2))/32.
!near bottom
j=j0 hf(j)=w21(2*j+1)*(12.*w(j0)-w(j0-1))/32.
!near topcycle_4: DO j=0,j0 VJP=w21(2*j+1)
m=1 IF(VJP>0.)m=-1 jp=j+m
IF(j==0.AND.m==-1 .OR. j==j0.AND.m==1)CYCLE cycle_4
VP=(w(j)+w(j+1))/2. dV=w(j+1)-w(j) dV1=w(jp+1)-w(ip)hf(j)=VJP*VP+ABS(VJP)*(-dV/2.+AMINMOD(dV1,4.*dV)/6.+AMINMOD(dV,4.*dV1)/3.) !hf=JVV_P
ENDDO cycle_4
FORALL(j=1:j0)Uq(2,i,j)=Uq(2,i,j)+hf(j)-hf(j-1)! +JVV_eta
ENDDO
ENDIF
! ***** Compute convection term at outlet using upwind-difference scheme
DO j=0,j0 i=if Uq(1,i,j)=0.
DO k=1,2
Uq(1,i,j) = Uq(1,i,j)+xix(k,2*i,2*j+1)&*(UJ(1,i,j)*(u(k,i,j)+u(k,i,j+1)-u(k,i-1,j)-u(k,i-1,j+1))/2. &+(UJX(2,i,j)+UJX(2,i,j+1))/2.*(u(k,i,j+1)-u(k,i,j)))ENDDO ENDDO
DO j=1,j0 i=i0 Uq(2,i,j)=0.
DO k=1,2
Uq(2,i,j) = Uq(2,i,j)+xix(k+2,2*i+1,2*j)&*((UJX(1,i,j)+UJX(1,i+1,j))/2.*(u(k,i+1,j)-u(k,i-1,j))/2.
&+UJ(2,i,j)*(u(k,i,j+1)+u(k,i+1,j+1)-u(k,i,j-1)-u(k,i+1,j-1))/4.)
ENDDO ENDDO! ***** Compute additional, pressure and diffusion terms
DO i=1,if
FORALL(j=0:jf)w2(j)=u(1,i,j)
CALL INTP(w2,w21,jf)!w21=u_U
FORALL(j=0:jf)w2(j)=u(2,i,j)
CALL INTP(w2,w22,jf)!w22=v_U
FORALL(j=0:jf)w2(j)=UJX(2,i,j) CALL INTP(w2,w23,jf)!w23=JV_U
FORALL(j=0:j0)w2(j)=(p(i-1,j)+p(i,j))/2. CALL DIFP(w2,w24,jf,1.)!w24=(p_eta)_U
DO j=0,j0 jj=2*j+1
ad = w21(jj)*(UJ(1,i,j)*dxix(i,j,1,1)+w23(jj)*dxix(i,j,1,2)) &!additional term+w22(jj)*(UJ(1,i,j)*dxix(i,j,1,3)+w23(jj)*dxix(i,j,1,4))
IF(i==if)ad = 0.
pr = gU(3,i,j)*(p(i,j)-p(i-1,j))-gU(2,i,j)*w24(j)!pressure termUq(1,i,j) = Uq(1,i,j)-ad+pr+(z(i,j+1)-z(i,j))/Re!residuals of NS eqn
rhs(1,i,j) = -dt*Uq(1,i,j)!rhs of linear eqns
ENDDO ENDDO
DO j=1,j0
FORALL(i=0:if)w1(i)=u(1,i,j)
CALL INTP(w1,w11,if)!w11=u_V41

42

非圧縮性流れの解法|MAC 型解法

FORALL(i=0:if)w1(i)=u(2,i,j)
CALL INTP(w1,w12,if)!w12=v_V
FORALL(i=0:if)w1(i)=UJX(1,i,j) CALL INTP(w1,w13,if)!w13=JU_V
FORALL(i=0:if)w1(i)=(p(i,j-1)+p(i,j))/2. CALL DIFP(w1,w14,if,1.)!w14=(p_xi)_V
i=i0 w14(i)=(w1(i+1)-w1(i-1))/2.
!outletDO i=0,i0 ii=2*i+1
ad = w11(ii)*(w13(ii)*dxix(i,j,2,1)+UJ(2,i,j)*dxix(i,j,2,2)) &!additional term+w12(ii)*(w13(ii)*dxix(i,j,2,3)+UJ(2,i,j)*dxix(i,j,2,4))IF(i==i0)ad = 0.
pr = -gV(2,i,j)*w14(i)+gV(1,i,j)*(p(i,j)-p(i,j-1))!pressure termUq(2,i,j) = Uq(2,i,j)-ad+pr-(z(i+1,j)-z(i,j))/Re!residuals of NS eqn
rhs(2,i,j) = -dt*Uq(2,i,j)!rhs of linear eqns
ENDDO ENDDO
IF(unsteady)THEN!only unsteadyIF(ma==1)FORALL(l=1:2,i=0:if,j=0:jf) UJ0(l,i,j)= UJ(l,i,j)
IF(ma==1)FORALL(l=1:2,i=0:if,j=0:jf)rhs0(l,i,j)=rhs(l,i,j)
FORALL(l=1:2,i=0:if,j=0:jf)rhs(l,i,j)=-(UJ(l,i,j)-UJ0(l,i,j))+(rhs0(l,i,j)+rhs(l,i,j))/2.
ENDIF! Compute UJ^* by implicit SMAC-scheme
tt=dt*theta
DO i=1,if DO j=0,j0
UJp=(UJ(1,i,j)+ABS(UJ(1,i,j)))/2. UJm=(UJ(1,i,j)-ABS(UJ(1,i,j)))/2.
VJU=(UJX(2,i,j)+UJX(2,i,j+1))/2.
VJp=(VJU+ABS(VJU))/2. VJm=(VJU-ABS(VJU))/2.
c(1,i,j)=-tt*(UJp+gU(3,i,j)/Re)!coefs of linear eqs
c(2,i,j)= tt*(UJm-gU(3,i,j)/Re) IF(i==if)c(2,i,j)=0.
c(3,i,j)=-tt*(VJp+gU(1,i,j)/Re)c(4,i,j)= tt*(VJm-gU(1,i,j)/Re)c(0,i,j)=aJ(2*i,2*j+1)-c(1,i,j)-c(2,i,j)-c(3,i,j)-c(4,i,j)
ENDDO ENDDO!Compute dU^* by modifies AF-method
DO j=0,j0 l=j+1 DO i=1,ifa(l,i,1)=c(1,i,j)a(l,i,2)=c(0,i,j)a(l,i,3)=c(2,i,j) b(l,i)=rhs(1,i,j)
ENDDO ENDDOCALL GAUSSZ(a,b,140,jf,if)
FORALL(j=0:j0,i=1:if)rhs(1,i,j)=b(j+1,i)!dU^**
DO j=0,j0 l=j+1 DO i=1,ifa(i,l,1)=c(3,i,j)a(i,l,2)=c(0,i,j)a(i,l,3)=c(4,i,j) b(i,l)=c(0,i,j)*rhs(1,i,j)
ENDDO ENDDOCALL GAUSSZ(a,b,140,if,jf)
FORALL(j=0:j0,i=1:if)rhs(1,i,j)=b(i,j+1)!dU^*! Compute VJ^* by implicit SMAC-scheme
DO i=0,i0 DO j=1,j0
UJV=(UJX(1,i,j)+UJX(1,i+1,j))/2.
UJp=(UJV+ABS(UJV))/2. UJm=(UJV-ABS(UJV))/2.
VJp=(UJ(2,i,j)+ABS(UJ(2,i,j)))/2. VJm=(UJ(2,i,j)-ABS(UJ(2,i,j)))/2.
c(1,i,j)=-tt*(UJp+gV(3,i,j)/Re)!coefs of linear eqs
c(2,i,j)= tt*(UJm-gV(3,i,j)/Re)c(3,i,j)=-tt*(VJp+gV(1,i,j)/Re)c(4,i,j)= tt*(VJm-gV(1,i,j)/Re)c(0,i,j)=aJ(2*i+1,2*j)-c(1,i,j)-c(2,i,j)-c(3,i,j)-c(4,i,j)
ENDDO ENDDO!Compute dV^* by modifies AF-method
DO j=1,j0 DO i=0,i0 k=i+1

バックステップ流路流れのプログラムa(j,k,1)=c(1,i,j)a(j,k,2)=c(0,i,j)a(j,k,3)=c(2,i,j) b(j,k)=rhs(2,i,j)
ENDDO ENDDOCALL GAUSSZ(a,b,140,j0,if)
FORALL(j=1:j0,i=0:i0)rhs(2,i,j)=b(j,i+1)!dV^**
DO j=1,j0 DO i=0,i0 k=i+1
a(k,j,1)=c(3,i,j)a(k,j,2)=c(0,i,j)a(k,j,3)=c(4,i,j) b(k,j)=c(0,i,j)*rhs(2,i,j)
ENDDO ENDDOCALL GAUSSZ(a,b,140,if,j0)
FORALL(j=1:j0,i=0:i0)rhs(2,i,j)=b(i+1,j)!dV^*
FORALL(j=0:j0,i=1:if)UJ(1,i,j)=UJ(1,i,j)+beta*aJ(2*i,2*j+1)*rhs(1,i,j) !JU^*FORALL(j=1:j0,i=0:i0)UJ(2,i,j)=UJ(2,i,j)+beta*aJ(2*i+1,2*j)*rhs(2,i,j) !JV^*FORALL(i=0:if,j=0:jf)f(i,j)=Uq(1,i,j)
lo=MAXLOC(f) FORALL(k=1:2)locf(k,1)=lo(k) fmax(1)=MAXVAL(f)lo=MINLOC(f) FORALL(k=1:2)locf(k,2)=lo(k) fmax(2)=MINVAL(f)FORALL(i=0:if,j=0:jf)f(i,j)=Uq(2,i,j)
lo=MAXLOC(f) FORALL(k=1:2)locf(k,3)=lo(k) fmax(3)=MAXVAL(f)lo=MINLOC(f) FORALL(k=1:2)locf(k,4)=lo(k) fmax(4)=MINVAL(f)ENDSUBROUTINE COMPU1! ********** Set up coefficients of difference equations for phi
SUBROUTINE COEF(co)PARAMETER(if=140,jf=25,if1=280,jf1=50)
DIMENSION co(0:if,0:jf,-1:1,-1:1),aJ(0:if1,0:jf1),gU(3,0:if,0:jf),gV(3,0:if,0:jf), &UJX(2,0:if,0:jf),w1(0:if,0:jf),gw(4,0:if,0:jf)
COMMON //na,Re,dt,i0,j0,lowRe,steady,unsteady
COMMON /METRIC2/aJ /METRIC3/gU,gV /COMPUZ1/UJX! ***** Set up coefficients of difference-equations for phi
DO i=0,if DO j=0,jf
gw(1,i,j)= dt*gU(3,i,j)gw(2,i,j)=-dt*gU(2,i,j)/4.
gw(3,i,j)=-dt*gV(2,i,j)/4.
gw(4,i,j)= dt*gV(1,i,j)
ENDDO ENDDOFORALL(l=-1:1,m=-1:1,i=0:if,j=0:jf)co(i,j,l,m)=0.
DO i=0,i0 DO j=0,j0!1st stepco(i,j,1, 1)= gw(2,i+1,j)co(i,j,0, 1)= gw(2,i+1,j)co(i,j,1, 0)= gw(1,i+1,j)co(i,j,0, 0)=-gw(1,i+1,j)co(i,j,1,-1)=-gw(2,i+1,j)co(i,j,0,-1)=-gw(2,i+1,j)
ENDDO
co(i, 0,1, 1)=co(i, 0,1, 1)+3.*gw(2,i+1, 0)!bottomco(i, 0,0, 1)=co(i, 0,0, 1)+3.*gw(2,i+1, 0)co(i, 0,1, 0)=co(i, 0,1, 0)-3.*gw(2,i+1, 0)co(i, 0,0, 0)=co(i, 0,0, 0)-3.*gw(2,i+1, 0)co(i,j0,1, 0)=co(i,j0,1, 0)+3.*gw(2,i+1,j0)!top
co(i,j0,0, 0)=co(i,j0,0, 0)+3.*gw(2,i+1,j0)
co(i,j0,1,-1)=co(i,j0,1,-1)-3.*gw(2,i+1,j0)
co(i,j0,0,-1)=co(i,j0,0,-1)-3.*gw(2,i+1,j0)
ENDDO
DO i=1,i0 DO j=0,j0!2nd stepco(i,j, 0, 1)=co(i,j, 0, 1)-gw(2,i,j)

43

44

非圧縮性流れの解法|MAC 型解法

co(i,j,-1, 1)=co(i,j,-1, 1)-gw(2,i,j)co(i,j, 0, 0)=co(i,j, 0, 0)-gw(1,i,j)co(i,j,-1, 0)=co(i,j,-1, 0)+gw(1,i,j)co(i,j, 0,-1)=co(i,j, 0,-1)+gw(2,i,j)co(i,j,-1,-1)=co(i,j,-1,-1)+gw(2,i,j)
ENDDO
co(i, 0, 0, 1)=co(i, 0, 0, 1)-3.*gw(2,i, 0)co(i, 0,-1, 1)=co(i, 0,-1, 1)-3.*gw(2,i, 0)co(i, 0, 0, 0)=co(i, 0, 0, 0)+3.*gw(2,i, 0)co(i, 0,-1, 0)=co(i, 0,-1, 0)+3.*gw(2,i, 0)co(i,j0, 0, 0)=co(i,j0, 0, 0)-3.*gw(2,i,j0)co(i,j0,-1, 0)=co(i,j0,-1, 0)-3.*gw(2,i,j0)co(i,j0, 0,-1)=co(i,j0, 0,-1)+3.*gw(2,i,j0)co(i,j0,-1,-1)=co(i,j0,-1,-1)+3.*gw(2,i,j0)
ENDDO
DO j=0,j0-1 DO i=0,i0
co(i,j, 1,1)=co(i,j, 1,1)+gw(3,i,j+1)
co(i,j, 1,0)=co(i,j, 1,0)+gw(3,i,j+1)
co(i,j, 0,1)=co(i,j, 0,1)+gw(4,i,j+1)
co(i,j, 0,0)=co(i,j, 0,0)-gw(4,i,j+1)
co(i,j,-1,1)=co(i,j,-1,1)-gw(3,i,j+1)
co(i,j,-1,0)=co(i,j,-1,0)-gw(3,i,j+1)
ENDDO
co(0,j,1,1)=co(0,j,1,1)+3.*gw(3,0,j+1)
co(0,j,1,0)=co(0,j,1,0)+3.*gw(3,0,j+1)
co(0,j,0,1)=co(0,j,0,1)-3.*gw(3,0,j+1)
co(0,j,0,0)=co(0,j,0,0)-3.*gw(3,0,j+1)
ENDDO
DO j=1,j0 DO i=0,i0
co(i,j, 1, 0)=co(i,j, 1, 0)-gw(3,i,j)co(i,j, 1,-1)=co(i,j, 1,-1)-gw(3,i,j)co(i,j, 0, 0)=co(i,j, 0, 0)-gw(4,i,j)co(i,j, 0,-1)=co(i,j, 0,-1)+gw(4,i,j)co(i,j,-1, 0)=co(i,j,-1, 0)+gw(3,i,j)co(i,j,-1,-1)=co(i,j,-1,-1)+gw(3,i,j)
ENDDO
co(0,j,1, 0)=co(0,j,1, 0)-3.*gw(3,0,j)co(0,j,1,-1)=co(0,j,1,-1)-3.*gw(3,0,j)co(0,j,0, 0)=co(0,j,0, 0)+3.*gw(3,0,j)co(0,j,0,-1)=co(0,j,0,-1)+3.*gw(3,0,j)
ENDDO
ENDSUBROUTINE COEF

!bottom!top!3rd step!inlet!4th step!inlet! ********** Compute static pressure p
SUBROUTINE COMPP(UJ,p,phi,co,resP,locf,fmax)
PARAMETER(if=140,jf=25,ife=70,jfe=13)!ife=(if+1)/2,jfe=(jf+1)/2
DIMENSION UJ(2,0:if,0:jf),p(0:if,0:jf),phi(0:if,0:jf),co(0:if,0:jf,-1:1,-1:1), &rhs(0:if,0:jf),a(ife,0:i0,3),b(ife,0:i0),r(ife),lo(2),locf(2,8),fmax(8)COMMON //na,Re,dt,i0,j0,lowRe,steady,unsteady
REAL(8) a,b
LOGICAL lowRe,steady,unsteady! Determine rhs of pressure difference eqn
FORALL(i=0:if,j=0:jf)phi(i,j)=0.
FORALL(i=0:i0,j=0:j0)rhs(i,j)=(UJ(1,i+1,j)-UJ(1,i,j)+UJ(2,i,j+1)-UJ(2,i,j)) !UJ=JU^*!IF(steady.AND.lowRe)THEN
IF(steady)THEN! ***** Solve pressure difference eqn by SOR method

バックステップ流路流れのプログラム

45

alp=1.5 bet=0.5!alp:over-relaxation, bet:damping
nf=4
n=0 100 n=n+1 resP=0.
ib=0 ie=i0 id=1 jb=0 je=j0 jd=1
!odd sweep, throughout region
IF(MOD(n,2)==0)THEN
ib=i0-1 ie=1 id=-1 jb=j0-1 je=1 jd=-1!even sweep, inside region
ENDIF
DO i=ib,ie,id DO j=jb,je,jd
res=-rhs(i,j)
DO ip=-1,1 ic=0 IF(i+ip==-1)ic=3
DO jp=-1,1 jc=0 IF(j+jp==-1)jc=3 IF(j+jp==if)jc=-3
res=res+co(i,j,ip,jp)*phi(i+ip+ic,j+jp+jc)!residuals
ENDDO ENDDO
phi(i,j)=phi(i,j)-alp*res/co(i,j,0,0)!over-relaxation of phi
resP=AMAX1(resP,ABS(res))!max residual
ENDDO ENDDO
IF(resP>1.E-5 .AND. n<nf)
GOTO 100!dicide convergence
ELSE
20 CONTINUE! ***** Solve pressure difference eqn by Tschebyscheff SLOR method
alp=1.25 bet=0.75
ifo=if-ife jfo=jf-jfe nf=4
n=0 110 n=n+1 resP=0.
! ***** Compute even i-columns
DO j=0,j0
DO ii=1,ife i=2*ii-2
a(ii,j,1)=co(i,j,0,-1) a(ii,j,2)=co(i,j,0, 0)*alp a(ii,j,3)=co(i,j,0, 1)b(ii,j) =rhs(i,j)-(1.-alp)*co(i,j,0,0)*phi(i,j)r(ii)=b(ii,j)ENDDO
DO jp=-1,1 jc=0 IF(j+jp==-1)jc=3 IF(j+jp==jf)jc=-3
DO ii=1,ife i=2*ii-2
im1=i-1 IF(im1==-1)im1=2
cop=co(i,j,-1,jp)*phi(im1,j+jp+jc)+co(i,j,1,jp)*phi(i+1,j+jp+jc)b(ii,j)=b(ii,j)-copr(ii)=r(ii)-cop-a(ii,j,jp+2)*phi(i,j+jp+jc)!residuals
ENDDO ENDDO
DO ii=1,ife
resP=AMAX1(resP,ABS(r(ii)))!max residual
ENDDO
ENDDO
CALL GAUSSP(a,b,ife,i0,j0,ife)
FORALL(j=0:j0,ii=1:ife)phi(2*ii-2,j)=b(ii,j)!even column phi! ***** Compute odd i-columns
DO j=0,j0
DO ii=1,ifo i=2*ii-1
a(ii,j,1)=co(i,j,0,-1) a(ii,j,2)=co(i,j,0, 0)*alp a(ii,j,3)=co(i,j,0, 1)b(ii,j) =rhs(i,j)-(1.-alp)*co(i,j,0,0)*phi(i,j)r(ii)=b(ii,j)ENDDO
DO jp=-1,1 jc=0 IF(j+jp==-1)jc=3 IF(j+jp==jf)jc=-3
DO ii=1,ifo i=2*ii-1
cop=co(i,j,-1,jp)*phi(i-1,j+jp+jc)+co(i,j,1,jp)*phi(i+1,j+jp+jc)b(ii,j)=b(ii,j)-copr(ii)=r(ii)-cop-a(ii,j,jp+2)*phi(i,j+jp+jc)!residuals
ENDDO ENDDO
DO ii=1,ifo

46

非圧縮性流れの解法|MAC 型解法

resP=AMAX1(resP,ABS(r(ii)))!max residual
ENDDO
ENDDO
CALL GAUSSP(a,b,ife,i0,j0,ifo)
FORALL(ii=1:ifo,j=0:j0)phi(2*ii-1,j)=b(ii,j)!odd column phi! Compute even j-rows
n=n+1 resP=0.
DO i=0,i0
DO jj=1,jfe j=2*jj-2
a(jj,i,1)=co(i,j,-1,0) a(jj,i,2)=co(i,j, 0,0)*alp a(jj,i,3)=co(i,j, 1,0)b(jj,i) =rhs(i,j)-(1.-alp)*co(i,j,0,0)*phi(i,j)r(jj)=b(jj,i)
ENDDO
DO ip=-1,1 ic=0 IF(i+ip==-1)ic=3
DO jj=1,jfe j=2*jj-2
jm1=j-1 IF(jm1==-1)jm1=2
jp1=j+1 IF(jp1==jf)jp1=j0-2
cop=co(i,j,ip,-1)*phi(i+ip+ic,jm1)+co(i,j,ip,1)*phi(i+ip+ic,jp1)b(jj,i)=b(jj,i)-copr(jj)=r(jj)-cop-a(jj,i,ip+2)*phi(i+ip+ic,j)!residuals
ENDDO ENDDO
DO jj=1,jfe
resP=AMAX1(resP,ABS(r(jj)))!max residual
ENDDO
ENDDO
CALL GAUSSP(a,b,ife,i0,i0,jfe)
FORALL(i=0:i0,jj=1:jfe)phi(i,2*jj-2)=b(jj,i)!even row phi! ***** Compute odd j-rows
DO i=0,i0
DO jj=1,jfo j=2*jj-1
a(jj,i,1)=co(i,j,-1,0) a(jj,i,2)=co(i,j, 0,0)*alp a(jj,i,3)=co(i,j, 1,0)b(jj,i) =rhs(i,j)-(1.-alp)*co(i,j,0,0)*phi(i,j)r(jj)=b(jj,i)
ENDDO
DO ip=-1,1 ic=0 IF(i+ip==-1)ic=3
DO jj=1,jfo j=2*jj-1
jm1=j-1 IF(jm1==-1)jm1=2
jp1=j+1 IF(jp1==jf)jp1=j0-2
cop=co(i,j,ip,-1)*phi(i+ip+ic,jm1)+co(i,j,ip,1)*phi(i+ip+ic,jp1)b(jj,i)=b(jj,i)-copr(jj)=r(jj)-cop-a(jj,i,ip+2)*phi(i+ip+ic,j)!residuals
ENDDO ENDDO
DO jj=1,jfo
resP=AMAX1(resP,ABS(r(jj)))!max residual
ENDDO
ENDDO
CALL GAUSSP(a,b,ife,i0,i0,jfo)
FORALL(jj=1:jfe,i=0:i0)phi(i,2*jj-1)=b(jj,i)!odd row phi
IF(resP>1.E-5 .AND. n<nf)
GOTO 110
ENDIF
FORALL(i=0:i0,j=0:j0)p(i,j)=p(i,j)+bet*phi(i,j)!p=p+bet*phi
lo=MAXLOC(phi) FORALL(k=1:2)locf(k,5)=lo(k) fmax(5)=MAXVAL(phi)lo=MINLOC(phi) FORALL(k=1:2)locf(k,6)=lo(k) fmax(6)=MINVAL(phi)ENDSUBROUTINE COMPP! ********** Compute volume flux JU and JVSUBROUTINE COMPU2(UJ,phi)バックステップ流路流れのプログラム
PARAMETER(if=140,jf=25)
DIMENSION UJ(2,0:if,0:jf),phi(0:if,0:jf),gU(3,0:if,0:jf),gV(3,0:if,0:jf), &pV(0:if),dpV(0:if),pU(0:jf),dpU(0:jf)
COMMON //na,Re,dt,i0,j0,lowRe,steady,unsteady
COMMON /METRIC3/gU,gV
FORALL(j=0:j0)phi(if,j)=-phi(i0,j)DO i=1,if!JU=JU^*-dt(g22phi_xi-g12phi_eta)
FORALL(j=0:j0)pU(j)=(phi(i-1,j)+phi(i,j))/2.
CALL DIFP(pU,dpU,jf,1.)
FORALL(j=0:j0) &UJ(1,i,j)=UJ(1,i,j)-dt*(gU(3,i,j)*(phi(i,j)-phi(i-1,j))-gU(2,i,j)*dpU(j))
ENDDO
DO j=1,j0!JV=JV^*-dt(-g12phi_xi+g11phi_eta)
FORALL(i=0:if)pV(i)=(phi(i,j-1)+phi(i,j))/2.
CALL DIFP(pV,dpV,if,1.)
i=i0 dpV(i)=(phi(i+1,j)-phi(i-1,j))/2.
!outlet
FORALL(i=0:i0) &UJ(2,i,j)=UJ(2,i,j)-dt*(-gV(2,i,j)*dpV(i)+gV(1,i,j)*(phi(i,j)-phi(i,j-1)))
ENDDO
ENDSUBROUTINE COMPU2! ********** Compute velocities u and v and vorticity zeta
SUBROUTINE COMPUZ(uJ,u,z,locf,fmax,CFLmax)
PARAMETER(if=140,jf=25,if1=280,jf1=50)
DIMENSION UJ(2,0:if,0:jf),u(2,0:if,0:jf),z(0:if,0:jf),div(0:if,0:jf),UJX(2,0:if,0:jf), &aJ(0:if1,0:jf1),gU(3,0:if,0:jf),gV(3,0:if,0:jf),xix(4,0:if1,0:jf1),w1(0:if), &w2(0:jf),w11(0:if1),w12(0:jf1),lo(2),locf(2,8),fmax(8)COMMON //na,Re,dt,i0,j0,lowRe,steady,unsteady
COMMON /COMPUZ1/UJX /METRIC2/aJ /METRIC3/gU,gV /METRIC4/xix! ***** Compute velocities at X
!u= xix4*JUX-xix2*JVX
DO i=0,if!v=-xix3*JUX+xix1*JVXFORALL(j=0:jf)w2(j)=UJ(1,i,j)
CALL INTUX(w2,w12,jf)
FORALL(j=0:jf)UJX(1,i,j)=w12(2*j)
ENDDO
DO j=0,jf
FORALL(i=0:if)w1(i)=UJ(2,i,j)
CALL INTVX(w1,w11,if)
FORALL(i=0:if)UJX(2,i,j)=w11(2*i)
ENDDO
FORALL(i=1:if,j=1:j0)u(1,i,j)= xix(4,2*i,2*j)*UJX(1,i,j)-xix(2,2*i,2*j)*UJX(2,i,j)
FORALL(i=1:if,j=1:j0)u(2,i,j)=-xix(3,2*i,2*j)*UJX(1,i,j)+xix(1,2*i,2*j)*UJX(2,i,j)! ***** Compute vorticity zeta at X
FORALL(i=0:if,j=0:jf)z(i,j)=0.
DO j=1,j0
FORALL(i=0:if)w1(i)=UJX(1,i,j)
CALL INTP(w1,w11,if)!w11=JU_V
FORALL(i=0:i0)w1(i)=gV(2,i,j)*w11(2*i+1)+gV(3,i,j)*UJ(2,i,j)
FORALL(i=1:i0)z(i,j)=w1(i)-w1(i-1)!*=(g21JU+g22JV)_,xi
i=0z(i,j)=-w1(i+2)+3.*w1(i+1)-2.*w1(i)!inlet boundi=if z(i,j)= w1(i-3)-3.*w1(i-2)+2.*w1(i-1)!outlet bound
ENDDO
DO i=0,if
FORALL(j=0:jf)w2(j)=UJX(2,i,j)
CALL INTVP(w2,w12,jf)!w12=JV_U
FORALL(j=0:j0)w2(j)=gU(1,i,j)*UJ(1,i,j)+gU(2,i,j)*w12(2*j+1)
FORALL(j=1:j0)z(i,j)=z(i,j)-(w2(j)-w2(j-1))!*=*-(g11JU+g12JV)_,eta

47

48

非圧縮性流れの解法|MAC 型解法

j=0
z(i,j)=-(-w2(j+1)/3.+3.*w2(j)) z(40,0)=0.
!bottom wall
j=jf z(i,j)=-( w2(j-2)/3.-3.*w2(j-1))!top wall
FORALL(j=0:jf)z(i,j)=z(i,j)/aJ(2*i,2*j)!zeta=*/J
ENDDO! ***** Divergence and max CFL number
CFLmax=0.
FORALL(i=0:i0,j=0:j0)div(i,j)=(UJ(1,i+1,j)-UJ(1,i,j) &+UJ(2,i,j+1)-UJ(2,i,j))/aJ(2*i+1,2*j+1)
lo=MAXLOC(div) FORALL(k=1:2)locf(k,7)=lo(k) fmax(7)=MAXVAL(div)lo=MINLOC(div) FORALL(k=1:2)locf(k,8)=lo(k) fmax(8)=MINVAL(div)
DO i=0,i0 DO j=0,j0
CFLmax=MAX(CFLmax,ABS(UJ(1,i,j)/aJ(2*i,2*j+1)),ABS(UJ(2,i,j)/aJ(2*i+1,2*j)))
ENDDO ENDDO
CFLmax=CFLmax*dt!max|CFL|ENDSUBROUTINE COMPUZ! ********** Compute location of reattachment point
SUBROUTINE REATTACH(x,UJ,UJX,if,jf,xrap,k)
DIMENSION x(2,0:if,0:jf),UJ(2,0:if,0:jf),UJX(2,0:if,0:jf),al(-3:8)
k=100 irap=45 10 irap=irap+1 IF(irap>105)RETURNIF(UJ(1,irap,0)<0.)
GOTO 10!irap: grid point near rap
ii=-4 12 ii=ii+1 i=irap+ii
du=UJ(1,i,0) d2u=UJX(1,i,1)-2.*du!du=DeltaJU, d2u=Delta^2JU
d3u=UJ(1,i,1)-3.*UJX(1,i,1)+3.*du!d3u=Delta^3JU
a=2. 13 fa = du+(a-1.)/2.*(d2u+(a-2.)/3.*d3u)!fa=f(alpha) by cubic interpolation
IF(ABS(fa)<.00001)
GOTO 14
dfa = d2u/2.+(2.*a-3.)/6.*d3u!dfa=f'(alpha)a = a-fa/dfa
GOTO 13!by Newton method
14 al(ii)=a IF(a>0.14)
GOTO 12!al(ii): alpha of u=0
IF(al(ii)<0.01)ii=ii-1!irap+ii: adjacent grid point to rap! using quadratic interpolation formula f(beta)=a0+da*beta+d2a*beta*beta/2.
! get beta from f(beta)=0 by Newton method
a2=al(ii-2) a1=al(ii-1) a0=al(ii)da=(3.*a0-4.*a1+a2)/2. d2a=a0-2.*a1+a2
k=0 b=0. 22 k=k+1 fb = a0+da*b+d2a*b*b/2.
!fb=f(beta)
IF(k>10)
GOTO 20
IF(ABS(fb)<0.01)
GOTO 21
dfb = da+d2a*b!dfb=f'(beta)b = b-fb/dfb
GOTO 22!by Newton method
21 xii=x(1,irap+ii,0) xi=x(1,irap+ii+1,0)20 xrap = xii+b*(xi-xii)!x of reattachment point
ENDSUBROUTINE REATTACH! ********** minmod limiter
FUNCTION AMINMOD(a,b)
s=SIGN(1.,a) AMINMOD=s*MAX(0.,MIN(ABS(a),s*b))
ENDFUNCTION! ********** Compute x- or y-differences at point P (pressure)
SUBROUTINE DIFP(u,du,n,h)
DIMENSION u(0:n),du(0:n)
FORALL(i=1:n-2)du(i)=(u(i+1)-u(i-1))/(2.*h)du(0)=(-3.*u(0)+4.*u(1)-u(2))/(2.*h)du(n-1)=(u(n-3)-4.*u(n-2)+3.*u(n-1))/(2.*h)
ENDSUBROUTINE DIFP
! ********** Compute x- or y-differences at point X (metrics)

バックステップ流路流れのプログラム
SUBROUTINE DIFX(u,du,n,h)
DIMENSION u(0:n),du(0:n)
FORALL(i=1:n-1)du(i)=(u(i+1)-u(i-1))/(2.*h)du(0)=(-3.*u(0)+4.*u(1)-u(2))/(2.*h)du(n)=(u(n-2)-4.*u(n-1)+3.*u(n))/(2.*h)
ENDSUBROUTINE DIFX
! ********** Solve in parallel, systems of linear eqns with modified tri-diagonal matrix! by Gaussian eliminaton
SUBROUTINE GAUSSP(a,b,ife,i0,n,if)
DIMENSION a(ife,0:i0,3),b(ife,0:i0)
REAL(8)a,b
FORALL(i=1:if)a(i,0,1)=a(i,0,1)/a(i,0,2)
FORALL(i=1:if)a(i,1,3)=a(i,1,3)-a(i,1,1)*a(i,0,1)
DO k=0,n-1 DO i=1,ifb(i,k)=b(i,k)/a(i,k,2) a(i,k,3)=a(i,k,3)/a(i,k,2)b(i,k+1)=b(i,k+1)-a(i,k+1,1)*b(i,k)
a(i,k+1,2)=a(i,k+1,2)-a(i,k+1,1)*a(i,k,3)
ENDDO ENDDOFORALL(i=1:if)a(i,n,1)=-a(i,n,3)*a(i,n-2,3)
FORALL(i=1:if)a(i,n,2)=a(i,n,2)-a(i,n,1)*a(i,n-1,3)
FORALL(i=1:if)b(i,n)=(b(i,n)-a(i,n,3)*b(i,n-2)-a(i,n,1)*b(i,n-1))/a(i,n,2)
DO k=n-1,0,-1 DO i=1,ifb(i,k)=b(i,k)-a(i,k,3)*b(i,k+1)
ENDDO ENDDO
DO i=1,if b(i,0)=b(i,0)-a(i,0,1)*b(i,2) ENDDO
ENDSUBROUTINE GAUSSP! ********** Solve in paralell, systems of linear eqns with tri-diagonal matrix! by Gaussian elimination
SUBROUTINE GAUSSZ(a,b,m1,m,n)
DIMENSION a(m1,m1,3),b(m1,m1)
DO k=1,n-1 DO l=1,mb(l,k) =b(l,k) /a(l,k,2)a(l,k,3)=a(l,k,3)/a(l,k,2)b(l,k+1) =b(l,k+1) -a(l,k+1,1)*b(l,k)
a(l,k+1,2)=a(l,k+1,2)-a(l,k+1,1)*a(l,k,3)
ENDDO ENDDOFORALL(l=1:m)b(l,n)=b(l,n)/a(l,n,2)
DO k=n-1,1,-1
FORALL(l=1:m)b(l,k)=b(l,k)-a(l,k,3)*b(l,k+1)
ENDDO
ENDSUBROUTINE GAUSSZ! ********** Interpolate u1 from u (U_P from U)
SUBROUTINE INTP(u,u1,n)DIMENSION u(0:n),u1(0:2*n)
FORALL(i=0:n)u1(2*i)=u(i)FORALL(i=1:n-2)u1(2*i+1)=(-u(i-1)+9.*(u(i)+u(i+1))-u(i+2))/16.
i=0u1(2*i+1)=(3.*u(i)+6.*u(i+1)-u(i+2))/8.
i=n-1 u1(2*i+1)=(3.*u(i+1)+6.*u(i)-u(i-1))/8.
ENDSUBROUTINE INTP
! ********** Interpolate U_X from U
SUBROUTINE INTUX(u,u1,jf)
DIMENSION u(0:jf),u1(0:2*jf)
FORALL(j=0:jf-1)u1(2*j+1)=u(j)

49

50

非圧縮性流れの解法|MAC 型解法

FORALL(j=2:jf-2)u1(2*j)=(-u(j-2)+9.*(u(j-1)+u(j))-u(j+1))/16.
j=0
u1(2*j)=0.
j=1
u1(2*j)=(3.*u(j-1)+2.*u(j)-u(j+1)/5.)/4.
j=jf-1 u1(2*j)=(3.*u(j)+2.*u(j-1)-u(j-2)/5.)/4.
j=jf
u1(2*j)=0.
ENDSUBROUTINE INTUX! ********** Interpolate V_P from V
SUBROUTINE INTVP(u,u1,n)DIMENSION u(0:n),u1(0:2*n)
FORALL(j=0:n)u1(2*j)=u(j)FORALL(j=1:n-2)u1(2*j+1)=(-u(j-1)+9.*(u(j)+u(j+1))-u(j+2))/16.
j=0
u1(2*j+1)=(12.*u(j+1)-u(j+2))/32.
j=n-1 u1(2*j+1)=(12.*u(j )-u(j-1))/32.
ENDSUBROUTINE INTVP! ********** Interpolate V_X from V
SUBROUTINE INTVX(u,u1,if)
DIMENSION u(0:if),u1(0:2*if)
FORALL(i=0:if-1)u1(2*i+1)=u(i)FORALL(i=2:if-2)u1(2*i)=(-u(i-2)+9.*(u(i-1)+u(i))-u(i+1))/16.
i=0u1(2*i)=0.
i=1
u1(2*i)=(3.*u(i-1)+2.*u(i)-u(i+1)/5.)/4.
i=if-1 u1(2*i)=(3.*u(i)+6.*u(i-1)-u(i-2))/8.
i=if
u1(2*i)=(3.*u(i-1)-u(i-2))/2.
ENDSUBROUTINE INTVX!********** Print out computational results in flow field
SUBROUTINE WRTF(f,z1,z2,mz,na)
PARAMETER(if=140,jf=25)
DIMENSION f(0:if,0:jf)
CHARACTER*4 z1,z2,form(7),area(7)DATA form/"(1H ","2X I","2,2X"," 36(","","F8.4",")) "/DATA area/" -3P"," -2P"," -1P",""," 1P"," 2P"," 3P"/
form(5)=area(mz+4)
OPEN(20,FILE='OUTPUT.dat')
WRITE(20,60)z1,z2,na
iif=(if-1)/35+1
DO ii=1,iif is=35*(ii-1) ie=MIN(35*ii,if)
DO j=jf,0,-1 WRITE(20,form)j,(f(i,j),i=is,ie) ENDDO
WRITE(20,61)(i,i=is,ie)
ENDDO
60 FORMAT(/1H 20X "* * * ", 2A4, " * * *", 10X "na =", I5, /)61 FORMAT(1H 4X 36I8//)ENDSUBROUTINE WRTF
! ********** Drawing of Grid and Computational Results
SUBROUTINE StepDF_CG(x,u,p,z,if,jf)
USE DFLIBDIMENSION x(2,0:if,0:jf),u(2,0:if,0:jf),p(0:if,0:jf),z(0:if,0:jf),pp(7:105,0:jf)
LOGICALstatusmode
TYPE(windowconfig) myscreen
TYPE(wxycoord)wxy
INTEGER(2) status,xwidth,yheight,fontnum,numfonts
INTEGER(4) oldcolor,colors(13)REAL(4)p0(0:if,0:jf),p1(0:if),p11(0:if),psi(0:if,0:jf),psik(13), &x1(0:70,0:25),y1(0:70,0:25),f1(0:70,0:25),&

バックステップ流路流れのプログラム51

x2(0:30,0:25),y2(0:30,0:25),f2(0:30,0:25),ak(13),Dps(13)REAL(8)xs,ys,xt,yt
CHARACTER*5 z1(7)myscreen.numxpixels = -1
myscreen.numypixels = -1
myscreen.numtextcols = -1
myscreen.numtextrows = -1
myscreen.numcolors= -1
myscreen.fontsize= -1
myscreen.title= " "Cstatusmode = SETWINDOWCONFIG(myscreen)
IF(.NOT. statusmode) statusmode = SETWINDOWCONFIG(myscreen)
statusmode =GETWINDOWCONFIG(myscreen)
xwidth = myscreen.numxpixels
yheight = myscreen.numypixels
oldcolor=SETBKCOLORRGB(#FFFFFF)CALL CLEARSCREEN($GCLEARSCREEN)! Set first window coordinates
CALL SETVIEWPORT(yheight/5.,INT2(0), yheight,yheight/5.)
status = SETWINDOW(.TRUE., -5.,-2., 11.,2.)! Draw partial computational grid
oldcolor=SETCOLORRGB(#606060)!gray
DO i=0,if xs=x(1,i,0) ys=x(2,i,0) CALL MOVETO_W(xs,ys,wxy)
DO j=1,jf xt=x(1,i,j) yt=x(2,i,j) status = LINETO_W(xt,yt) ENDDO
ENDDO
DO j=0,jf xs=x(1,0,j) ys=x(2,0,j) CALL MOVETO_W(xs,ys,wxy)
DO i=1,if xt=x(1,i,j) yt=x(2,i,j) status = LINETO_W(xt,yt) ENDDO
ENDDO
oldcolor=SETCOLORRGB(#000000)
CALL CHARAC(1.0,-1.8,'COMPUTATIONAL GRID')! Set second window coordinates
OPEN(20,FILE='OUTPUT.dat')
CALL SETVIEWPORT(yheight/5.,yheight/5., yheight,3.*yheight/5.)
status = SETWINDOW(.TRUE., -6.,-8., 18.,4.)! Set colors
colors(1) =#E070A0 colors(2) =#FF7070!bluecolors(3) =#D0A070 colors(4) =#A0D070 colors(5) =#70FF70 !green
colors(6) =#60FFA0 colors(7) =#50FFD0 colors(8) =#40FFFF !yellow
colors(9) =#50E0FF colors(10)=#60C0FF colors(11)=#7090FF
colors(12)=#7070FF!redcolors(13)=#A070E0
status = REMAPALLPALETTERGB(colors)! Draw pressure distribution
DO j=0,jf-1!interpolate p_X
FORALL(i=0:if-1)p1(i)=p(i,j)
CALL INTPX(p1,p11,if,if)
FORALL(i=0:if)p0(i,j)=p11(i)ENDDO
DO i=0,if
FORALL(j=0:jf-1)p1(j)=p0(i,j)
CALL INTPX(p1,p11,if,jf)
FORALL(j=0:jf)p0(i,j)=p11(j)ENDDODATA Dps/.005, .00667, .0075, .010, .0125, .015, .01667, .02, .025, &.03, .03333, .04, .05/
!standard pressure increments

52

非圧縮性流れの解法|MAC 型解法

FORALL(i=7:105,j=0:jf)pp(i,j)=p0(i,j)
pmax=MAXVAL(pp) pmin=MINVAL(pp)Dp=(pmax-pmin)/13.
IF(Dp<.004 .OR. Dp>.06)
RETURN!no suitable Dps
k=0 10 k=k+1 IF(Dps(k)<Dp .AND. k<13)
GOTO 10!determine Dps
kkk=k Dpk=Dps(k) imin=INT(pmin/Dpk) IF(pmin>0.)imin=imin+1!No of min pressureak(12)=Dpk*(FLOAT(imin)+11.5)
11 IF(ak(12)>=pmax)THEN imin=imin-1 ak(12)=ak(12)-Dpk GOTO 11!determine imin
ENDIF
FORALL(k=1:12)ak(k)=Dpk*(FLOAT(imin+k)-.5)
FORALL(i=0:70,j=0:jf)x1(i,j)=x(1,i+6,j) y1(i,j)=x(2,i+6,j) f1(i,j)=p0(i+6,j) ENDFORALLCALL DISTRIB(x1,y1,f1,70,jf,ak,colors,1,13)
WRITE(20,'(5X A6,F8.4,5X A6,F8.4)')'pmax =',pmax,'pmin =',pmin! Draw velocity vectors
oldcolor=SETCOLORRGB(#000000)!blackamag=.6
DO i=7,75,4 DO j=1,jf-1
CALL LINE(x(1,i,j),x(2,i,j),x(1,i,j)+amag*u(1,i,j),x(2,i,j)+amag*u(2,i,j))
ENDDO ENDDO
status = SETWINDOW(.TRUE., 18.,-3., 42.,9.)! Draw pressure distribution
FORALL(i=0:30,j=0:jf)x2(i,j)=x(1,i+75,j) y2(i,j)=x(2,i+75,j) f2(i,j)=p0(i+75,j) ENDFORALLCALL DISTRIB(x2,y2,f2,30,jf,ak,colors,1,13)! Draw velocity vectors
oldcolor=SETCOLORRGB(#000000)
DO i=75,99,4 DO j=1,jf-1
CALL LINE(x(1,i,j),x(2,i,j),x(1,i,j)+amag*u(1,i,j),x(2,i,j)+amag*u(2,i,j))
ENDDO ENDDO! Color-pressure relation
CALL COLORVAL(39.0,-0.95,39.3,2.95,colors,13,kkk,Dpk,imin)
oldcolor=SETCOLORRGB(#404040)
CALL CHARAC(23.0,-2.3,'PRESSURE DISTRIBUTION AND VELOCITY VECTORS')! Set third window coordinates
CALL SETVIEWPORT(yheight/5.,3.*yheight/5., yheight,yheight)
status = SETWINDOW(.TRUE., -6.,-8., 18.,4.)! Draw vorticity distribution
FORALL(i=0:70,j=0:jf)x1(i,j)=x(1,i+6,j) y1(i,j)=x(2,i+6,j) f1(i,j)=z(i+6,j) ENDFORALL
FORALL(k=1:13)ak(k)=-2.6+.4*FLOAT(k)!paint colors(k) between ak(k-1) and ak(k)CALL DISTRIB(x1,y1,f1,70,jf,ak,colors,1,13)! Draw streamlines
DO i=0,if
psi(i,0)=0.
DO j=1,jf
psi(i,j)=psi(i,j-1)+(u(1,i,j-1)+u(1,i,j))/2.*(x(2,i,j)-x(2,i,j-1)) &-(u(2,i,j-1)+u(2,i,j))/2.*(x(1,i,j)-x(1,i,j-1))
ENDDO
dpsi=3.-psi(i,jf)
FORALL(j=1:jf)psi(i,j)=psi(i,j)+dpsi*FLOAT(j)/jf
ENDDO
FORALL(i=0:70,j=0:jf)f1(i,j)=psi(i+6,j)oldcolor=SETCOLORRGB(#606060)!gray
FORALL(k=1:13)psik(k)=.25*FLOAT(k-1)
CALL CONTOUR(x1,y1,f1,70,jf,psik,13)

バックステップ流路流れのプログラム
oldcolor=SETCOLORRGB(#808080)!light gray
FORALL(k=1:13)psik(k)=-.01*FLOAT(k)CALL CONTOUR(x1,y1,f1,70,jf,psik,13)
oldcolor=SETCOLORRGB(#606060)!gray
CALL LINE(-6.0,0.0,0.0,0.0)
CALL LINE(0.0,0.0,0.0,-1.0)
CALL LINE(0.0,-1.0,18.0,-1.0) CALL LINE(-6.0,3.0,18.0,3.0)
status = SETWINDOW(.TRUE., 18.,-3.1, 42.,9)! Draw vorticity distribution
FORALL(i=0:30,j=0:jf)x2(i,j)=x(1,i+75,j) y2(i,j)=x(2,i+75,j) f2(i,j)=z(i+75,j) ENDFORALLCALL DISTRIB(x2,y2,f2,30,jf,ak,colors,1,13)! Draw streamlines
FORALL(i=0:30,j=0:jf)f2(i,j)=psi(i+75,j)oldcolor=SETCOLORRGB(#606060)!gray
FORALL(k=1:13)psik(k)=.25*FLOAT(k-1)
CALL CONTOUR(x2,y2,f2,30,jf,psik,13)
oldcolor=SETCOLORRGB(#808080)!light gray
FORALL(k=1:13)psik(k)=-.01*FLOAT(k)CALL CONTOUR(x2,y2,f2,30,jf,psik,13)
oldcolor=SETCOLORRGB(#606060)!gray
CALL LINE(18.0,-1.0,42.0,-1.0) CALL LINE(18.0,3.0,42.0,3.0)! Color-vorticity relation
xs=39.0 xt=39.3
DO k=1,13
oldcolor=SETCOLORRGB(colors(k))
ys=-0.95+0.3*FLOAT(k) yt=-0.95+0.3*FLOAT(k-1)
status=RECTANGLE_W($GFILLINTERIOR,xs,ys,xt,yt)
ENDDO
oldcolor=SETCOLORRGB(#000000)
ys=2.95 yt=-.95 status = RECTANGLE_W($GBORDER,xs,ys,xt,yt)
numfonts=INITIALIZEFONTS()
fontnum =SETFONT('t''Arial''h14w6')
DATA z1/'-2.4 ','-1.6 ','-0.8 ',' 0.0 ',' 0.8 ',' 1.6 ',' 2.4 '/DO k=1,7 xt=39.4 yt=-0.55+0.6*FLOAT(k-1)
CALL MOVETO_W(xt,yt,wxy) CALL OUTGTEXT(z1(k))
ENDDOCALL CHARAC(24.0,-2.3,'VORTICITY DISTRIBUTION AND STREAMLINES')ENDSUBROUTINE StepDF_CG! ********** character
SUBROUTINE CHARAC(x,y,A)
USE DFLIBTYPE (windowconfig) myscreen
TYPE (wxycoord)wxy
INTEGER(2)numfonts,fontnum
REAL(8)xt,yt
CHARACTER*100
A
numfonts=INITIALIZEFONTS()
fontnum=SETFONT("t'Arial''h16w6'")
xt=x yt=y+0.33 CALL MOVETO_W(xt,yt,wxy)
ENDSUBROUTINE CHARAC

CALL OUTGTEXT(A)! ********** Values of colors
SUBROUTINE COLORVAL(x1,y1,x2,y2,colors,kf,kkk,Dpk,imin)
USE DFLIB
TYPE(wxycoord) wxy
INTEGER*2
status,numfonts,fontnum,intv(13)

53

54

非圧縮性流れの解法|MAC 型解法

INTEGER*4
oldcolor,colors(kf)REAL*8xs,ys,xt,yt
CHARACTER*5 z1(-20:50)
dy=(y2-y1)/FLOAT(kf)xs=x1 xt=x2
DO k=1,kf
oldcolor=SETCOLORRGB(colors(k))
ys=y1+dy*FLOAT(k) yt=y1+dy*FLOAT(k-1)
status = RECTANGLE_W($GFILLINTERIOR,xs,ys,xt,yt)
ENDDO
oldcolor=SETCOLORRGB(#000000)
ys=y2 yt=y1 status = RECTANGLE_W($GBORDER,xs,ys,xt,yt)
numfonts=INITIALIZEFONTS()
fontnum =SETFONT('t''Arial''h14w6')
DATA z1/'-0.20','-0.19','-0.18','-0.17','-0.16','-0.15','-0.14','-0.13','-0.12','-0.11', &'-0.10','-0.09','-0.08','-0.07','-0.06','-0.05','-0.04','-0.03','-0.02','-0.01', &' 0.00',' 0.01',' 0.02',' 0.03',' 0.04',' 0.05',' 0.06',' 0.07',' 0.08',' 0.09', &' 0.10',' 0.11',' 0.12',' 0.13',' 0.14',' 0.15',' 0.16',' 0.17',' 0.18',' 0.19', &' 0.20',' 0.21',' 0.22',' 0.23',' 0.24',' 0.25',' 0.26',' 0.27',' 0.28',' 0.29', &' 0.30',' 0.31',' 0.32',' 0.33',' 0.34',' 0.35',' 0.36',' 0.37',' 0.38',' 0.39', &' 0.40',' 0.41',' 0.42',' 0.43',' 0.44',' 0.45',' 0.46',' 0.47',' 0.48',' 0.49', &' 0.50'/!pressure value of each color band
DATA intv/2, 3, 4, 2, 4, 2, 3, 3*2, 3, 2*2/
!text intervals
intvk=intv(kkk)IF(imin>=0)THEN itextmin=((imin+intvk-1)/intvk)*intvk
ELSE
itextmin=( imin
/intvk)*intvk ENDIF!No of bottom text
l=0 12 l=l+1
xt=x2+0.1 yt=y1+dy*(itextmin-imin+intvk*(l-1))+0.4!coordinates of character
itext=INT((itextmin+intvk*(l-1))*(Dpk+.0001)*100.)!No of press character
CALL MOVETO_W(xt,yt,wxy) CALL OUTGTEXT(z1(itext))IF(itextmin-imin+intvk*l<=12)
GOTO 12ENDSUBROUTINE COLORVAL! ********** Draw contours of f with oldcolor
SUBROUTINE CONTOUR(x,y,f0,if,jf,fk,kf)
USE DFLIBDIMENSION x(0:if,0:jf),y(0:if,0:jf),f0(0:if,0:jf),fk(kf)
TYPE(wxycoord) wxy
INTEGER(2)statusREAL(4)f(0:if,0:jf)
REAL(8)xs,ys
e=.00001
DO k=1,kf
DO i=0,if DO j=0,jf
f(i,j)=f0(i,j)-fk(k)IF(ABS(f(i,j))<e) f(i,j)=SIGN(e,f(i,j))
ENDDO ENDDO
DO 10 i=0,if-1 DO 10 j=0,jf-1 imoveto=0!for all elements
i1=i
j1=j
x1=x(i1,j1) y1=y(i1,j1) f1=f(i1,j1)
i2=i+1 j2=jx2=x(i2,j2) y2=y(i2,j2) f2=f(i2,j2)
i3=i
j3=j+1 x3=x(i3,j3) y3=y(i3,j3) f3=f(i3,j3)
i4=i+1 j4=j+1 x4=x(i4,j4) y4=y(i4,j4) f4=f(i4,j4)
IF(f1*f2<0)THEN
xs=(f1*x2-f2*x1)/(f1-f2) ys=(f1*y2-f2*y1)/(f1-f2)!starting point
CALL MOVETO_W(xs,ys,wxy) imoveto=1
ENDIF
IF(f3*f4<0)THENバックステップ流路流れのプログラム
xs=(f3*x4-f4*x3)/(f3-f4) ys=(f3*y4-f4*y3)/(f3-f4)
IF(imoveto==1)THEN
status = LINETO_W(xs,ys)
imoveto=0
ELSE CALL MOVETO_W(xs,ys,wxy) imoveto=1
ENDIF
ENDIF
IF(f1*f3<0)THEN
xs=(f1*x3-f3*x1)/(f1-f3) ys=(f1*y3-f3*y1)/(f1-f3)
IF(imoveto==1)THEN
status = LINETO_W(xs,ys)
imoveto=0
GOTO 10ELSE CALL MOVETO_W(xs,ys,wxy) imoveto=1
ENDIF
ENDIF
IF(imoveto==1)THEN
xs=(f2*x4-f4*x2)/(f2-f4) ys=(f2*y4-f4*y2)/(f2-f4)
status = LINETO_W(xs,ys)
ENDIF
10 CONTINUEENDDOENDSUBROUTINE CONTOUR! ********** Draw distribution of f with colors
SUBROUTINE DISTRIB(x,y,f,if,jf,ak,colors,kb,kf)
USE DFLIBDIMENSION x(0:if,0:jf),y(0:if,0:jf),f(0:if,0:jf),ak(kb:kf),colors(kb:kf)
TYPE(wxycoord)
wxy, poly(6)INTEGER(2) status,m(4)INTEGER(4) oldcolor,colors
REAL(4)fe(4)REAL(8)xs,ys,xt,yt,xs0,ys0,xt0,yt0,xs1,ys1,xt1,yt1,wx,wy
jp=jf+1 ak(kf)=1.E10
DO i=0,if-1 DO j=0,jf-1 m1=i*jp+j+1!for all quadrilateral elementsm(1)=m1
fe(1)=f(i,j)m(2)=m1+jp
fe(2)=f(i+1,j)m(3)=m1+1
fe(3)=f(i,j+1)m(4)=m1+jp+1 fe(4)=f(i+1,j+1)! m() and fe() are four node numbers and function values of each element, respectively.
! l1 l2 l3 or l4 is one of four element node numbers 1 2 3 4, and the fe-values increase in! order. First locate max and min fe-value nodes l1 and l4, then determine middle fe-value! nodes l2 and l3.
l1=MINLOC(fe,DIM=1) l4=MAXLOC(fe,DIM=1)
sum=FLOAT(10-l1-l4) prd=FLOAT(24/l1/l4) sqr=SQRT(sum*sum-4.*prd)
l2=INT((sum-sqr)/2.+.1) l3=INT((sum+sqr)/2.+.1)IF(fe(l2)>fe(l3))THEN l=l2 l2=l3 l3=l ENDIF! a1<=a2<=a3<=a4 are nodal values of function f.
i1=(m(l1)-1)/jp j1=MOD(m(l1)-1,jp) x1=x(i1,j1) y1=y(i1,j1) a1=f(i1,j1)
i2=(m(l2)-1)/jp j2=MOD(m(l2)-1,jp) x2=x(i2,j2) y2=y(i2,j2) a2=f(i2,j2)
i3=(m(l3)-1)/jp j3=MOD(m(l3)-1,jp) x3=x(i3,j3) y3=y(i3,j3) a3=f(i3,j3)
i4=(m(l4)-1)/jp j4=MOD(m(l4)-1,jp) x4=x(i4,j4) y4=y(i4,j4) a4=f(i4,j4)
k1=kb-1 11 k1=k1+1 IF(ak(k1)<a1)
GOTO 11k2=kb-1 12 k2=k2+1 IF(ak(k2)<a2)
GOTO 12k3=kb-1 13 k3=k3+1 IF(ak(k3)<a3)
GOTO 13k4=kb-1 14 k4=k4+1 IF(ak(k4)<a4)
GOTO 14! ak(k1)<=ak(k2)<=ak(k3)<=ak(k4) Small range ak(k-1)-ak(k) is painted by colors(k).
! i.e. -infty-ak(kb) by colors(kb), ak(kf-1)-ak(kf):infty by colors(kf).
! k1<=k2<=k3<=k4. ntype shows arrangement of k1 k2 k3 k4,!
ntype=1: k1-k2-k4-k3-k1

55

56
!
!

非圧縮性流れの解法|MAC 型解法

ntype=2: k1-k2-k3-k4-k1,
ntype=3: k1-k3-k2-k4-k1(saddle point).
ntype=1
LA = i3==i1.OR.j3==j1 IF(LA==.FALSE.)ntype=2
LB = i2==i1.OR.j2==j1 IF(LB==.FALSE.)ntype=3! xs ys and xt yt are coordinates of points on edge of quadrilateral element, whose function! value is a. The coordinates are determined by linear interpolation.
DO 20 k=k1,k4 a=ak(k)IF(ntype==3)THEN nvertx=4!temporarypoly(1)%wx=x1 poly(1)%wy=y1
poly(2)%wx=x3 poly(2)%wy=y3
poly(3)%wx=x2 poly(3)%wy=y2
poly(4)%wx=x4 poly(4)%wy=y4
GOTO 21
ENDIF
IF(k==k4)
GOTO 22
IF(ntype==1)THEN!ntype=1
IF(k<k3)THEN xs=(a1-a)/(a1-a3)*x3+(a-a3)/(a1-a3)*x1 !k1<=k<k3
ys=(a1-a)/(a1-a3)*y3+(a-a3)/(a1-a3)*y1
ELSE
xs=(a3-a)/(a3-a4)*x4+(a-a4)/(a3-a4)*x3 !k3<=k<k4
ys=(a3-a)/(a3-a4)*y4+(a-a4)/(a3-a4)*y3
ENDIF
IF(k<k2)THEN xt=(a1-a)/(a1-a2)*x2+(a-a2)/(a1-a2)*x1 !k1<=k<k2
yt=(a1-a)/(a1-a2)*y2+(a-a2)/(a1-a2)*y1
ELSE
xt=(a2-a)/(a2-a4)*x4+(a-a4)/(a2-a4)*x2 !k2<=k<k4
yt=(a2-a)/(a2-a4)*y4+(a-a4)/(a2-a4)*y2
ENDIF
ELSEIF(ntype==2)THEN!ntype=2
xs=(a1-a)/(a1-a4)*x4+(a-a4)/(a1-a4)*x1 !k1<=k<k4
ys=(a1-a)/(a1-a4)*y4+(a-a4)/(a1-a4)*y1
IF(k<k2)THEN xt=(a1-a)/(a1-a2)*x2+(a-a2)/(a1-a2)*x1 !k1<=k<k2
yt=(a1-a)/(a1-a2)*y2+(a-a2)/(a1-a2)*y1
ELSEIF(k<k3)THEN
xt=(a2-a)/(a2-a3)*x3+(a-a3)/(a2-a3)*x2 !k2<=k<k3
yt=(a2-a)/(a2-a3)*y3+(a-a3)/(a2-a3)*y2
ELSE
xt=(a3-a)/(a3-a4)*x4+(a-a4)/(a3-a4)*x3 !k3<=k<k4
yt=(a3-a)/(a3-a4)*y4+(a-a4)/(a3-a4)*y3
ENDIF
ELSE
xs=(a1-a)/(a1-a4)*x4+(a-a4)/(a1-a4)*x1 !ntype=3
ys=(a1-a)/(a1-a4)*y4+(a-a4)/(a1-a4)*y1
IF(k2<k3.AND.k2<=k.AND.k<=k3)THEN
xs1=(a2-a)/(a2-a3)*x3+(a-a3)/(a2-a3)*x2
ys1=(a2-a)/(a2-a3)*y3+(a-a3)/(a2-a3)*y2
ENDIF
IF(k<k3)THEN xt=(a1-a)/(a1-a3)*x3+(a-a3)/(a1-a3)*x1
yt=(a1-a)/(a1-a3)*y3+(a-a3)/(a1-a3)*y1
ENDIF
IF(k2<=k)THEN xt1=(a2-a)/(a2-a4)*x4+(a-a4)/(a2-a4)*x2
yt1=(a2-a)/(a2-a4)*y4+(a-a4)/(a2-a4)*y2
ENDIF
ENDIF
22 CONTINUE! Divide into some polygons quadrilateral element in accordance with stepped function a(k).
! poly() are coordinates of vertices of polygon, and nvertx is number of the vertices.
! Coloring is performed in order k1...k4 in each element, but the following program written! in order k1 k4 k1<k<k4.
IF(k==k1)THEN!k==k1poly(1)%wx=x1 poly(1)%wy=y1

バックステップ流路流れのプログラム
poly(2)%wx=xs poly(2)%wy=yspoly(3)%wx=xt poly(3)%wy=yt nvertx=3
IF(k==k2)THEN
poly(4)%wx=x2 poly(4)%wy=y2 nvertx=4
ENDIF
IF(k==k4)THEN
IF(ntype==1)THEN
poly(2)%wx=x3 poly(2)%wy=y3
poly(3)%wx=x4 poly(3)%wy=y4
ENDIF
IF(ntype==2)THEN
poly(2)%wx=x4 poly(2)%wy=y4
poly(3)%wx=x3 poly(3)%wy=y3
ENDIF
ELSEIF(k==k3)THEN
poly(5)%wx=x2 poly(5)%wy=y2 nvertx=5
IF(ntype==1)THEN
poly(2)%wx=x3 poly(2)%wy=y3
poly(3)%wx=xs poly(3)%wy=yspoly(4)%wx=xt poly(4)%wy=yt
ENDIF
IF(ntype==2)THEN
poly(4)%wx=x3 poly(4)%wy=y3
ENDIF ENDIF
ELSEIF(k==k4)THEN
poly(1)%wx=xt0 poly(1)%wy=yt0
poly(2)%wx=xs0 poly(2)%wy=ys0
poly(3)%wx=x4
poly(3)%wy=y4 nvertx=3
IF(k==k3)THEN
nvertx=4
IF(ntype==1)THEN
poly(3)%wx=x3 poly(3)%wy=y3
poly(4)%wx=x4 poly(4)%wy=y4
ENDIF
IF(ntype==2)THEN
poly(4)%wx=x3 poly(4)%wy=y3
ENDIF ENDIF
IF(k==k2)THEN
poly(5)%wx=x2 poly(5)%wy=y2 nvertx=5
ENDIF
ELSE
poly(1)%wx=xt0 poly(1)%wy=yt0
poly(2)%wx=xs0 poly(2)%wy=ys0
poly(3)%wx=xs
poly(3)%wy=yspoly(4)%wx=xt
poly(4)%wy=yt nvertx=4
IF(k==k2.AND.k<k3)THEN
poly(5)%wx=x2 poly(5)%wy=y2 nvertx=5
ENDIF
IF(k==k3)THEN
nvertx=5
IF(ntype==1)THEN
poly(3)%wx=x3 poly(3)%wy=y3
poly(4)%wx=xs poly(4)%wy=yspoly(5)%wx=xt poly(5)%wy=yt
ENDIF
IF(ntype==2)THEN
poly(5)%wx=x3 poly(5)%wy=y3
ENDIF
IF(k==k2)THEN!k==k1<k2!k==k1==k2<k3!k==k1==k2==k3==k4!k==k1==k2==k3<k4!k1<k==k4!k3<k==k4!k2<k==k3==k4!k1<k==k2==k3==k4!k1<k<k4!k/=k2,k3!k1<k==k2<k3!k2<k==k3<k4!k1<k==k2==k3<k4

57

58

非圧縮性流れの解法|MAC 型解法

poly(6)%wx=x2 poly(6)%wy=y2 nvertx=6
ENDIF ENDIF
ENDIF
IF(k<k4)THEN
xs0=xs ys0=ys xt0=xt yt0=yt
ENDIF
21 oldcolor=SETCOLORRGB(colors(k))
status=POLYGON_W($GFILLINTERIOR,poly,nvertx)
20 CONTINUEENDDO ENDDO
ENDSUBROUTINE DISTRIB! ********** Interpolate values at U from values at P
SUBROUTINE INTPX(u,u1,if,n)
DIMENSION u(0:if),u1(0:if)
FORALL(i=2:n-2)u1(i)=(-u(i-2)+9.*(u(i-1)+u(i))-u(i+1))/16.
i=0  u1(i)=(3.*u(i)-u(i+1))/2.
i=1  u1(i)=(3.*u(i-1)+6.*u(i)-u(i+1))/8.
i=n-1 u1(i)=(3.*u(i)+6.*u(i-1)-u(i-2))/8.
i=n  u1(i)=(3.*u(i-1)-u(i-2))/2.
ENDSUBROUTINE INTPX! ********** Draw line
SUBROUTINE LINE(x1,y1,x2,y2)
USE DFLIBTYPE (wxycoord) wxy
INTEGER(2)statusREAL(8)xs,ys,xt,yt
xs=x1 ys=y1 xt=x2 yt=y2 CALL MOVETO_W(xs,ys,wxy) status = LINETO_W(xt,yt)
ENDSUBROUTINE LINE

次にこのプログラムのおもなサブルーチンを説明する．サブルーチン GRID は，楕円型微分方程式の境界値問題を解いて曲線座標格子を形成する解析的格子形成法 (analytical grid generation) のプログラムである．計算はすべて写像面上の単位正方形格子に対して行われる．バックステップ 流路は写像面上の長方形
領域 025 に写像される．格子間隔  =  = 1．上流境界 (x = ;10) の格子点番号140 0 

は i = 0，下流境界 (x = 70) は if = 140，ステップ (x = 0) の角は i1 = 35，隅は i2 = 40 である．また下方壁面 (y = 0 (xjf = 25 である．

0) x = 0 (;1 < y < 0) y = ;1 (x  0)) の格子点番号は j = 0，上方壁面 (y = 3) は曲線座標格子は基本的に次の Poisson 方程式の境界値問題を解くことによって形成される．

xx + yy = P

xx + yy = Q(13.84)

ただし P Q は制御関数と呼ばれ主に格子間隔の制御に使われる．しかしながらこの方程式は，独立変数がxy で  面上での計算に適しない．格子形成の計算には上式を書換えた次式が用いられる．
L(x) = ;J (x P + x Q)
L(y) = ;J (y P + y Q)

ただし L は微分演算子で
2
2@2
L = g~22 @@ 2 ; 2~g12 @ @@ +~g11 @
2(13.85a)(13.85b)59

バックステップ流路流れのプログラム

g~ij = gij =J ，ヤコビアン J と測度テンソル gij は前節で定義されたものである19 ．

式 (13.85) の境界値問題の境界条件は次のように与えられる．下方境界には境界点の座標が等比級数に;10，下流境界点の x = 70 が与えられる．上方境
界点の x の値は適当に，また内点の xy の値はとりあえず x ij = (1 ; j ) x i0 + j x ijf によって補間される．
 ) = ;~3 +1:6 ~2 +0:4 ~ ~ = =f である．なおこの 3 次式は条件 (0) = 0 (1) = 1  (0)ただし (~
よって，また上方境界点の y = 3，上流境界点の x =0= 0:4  (1) = 0:6 を満足するように作られたものである．
0

制御関数 P は，ここでは下方境界上に与えられた格子間隔が領域内部まで適当に浸透するように次式によって与えられる．

Pij = Pi0 or = t;  Pi0t
Pi0 = ( xx )i0 = x ;2 x
i+10

1

xi+10i 10;1;x ; x ;x
i0

i0



i 10
;

また制御関数 Q は，j 方向の格子間隔が流路中央で粗く境界とりわけ下方境界で細かくなるように次式によって与えられる．

dy~ 2 ~ (~y) d
Qij = (yy )ij =  (y)ij = dy
j d~
00

00

25 
dy~ 2 d =dy ij d~ (yijf ; yi0 )2ij

~ (~y) = 8~y ; 500

ただし 下の行の第 1 式は ~ = =25 y~ = (y; y0)=(yjf ; y0) なる関係から求められる倍率，また第 2 式は3 次式 ~(~y) = (8~y ; 15~y +13~y)=6 (0 y~ 1) の 2 階微分である．この３次式は条件 ~(0) = 0 ~(1) =1 ~ (0) = 13=6 ~ (1) = 7=6 を満足するもので， 方向の格子間隔の分布はこの式によって規定される．な
3

0

2

0

おプログラムには，格子間隔の差の小さいものから大きいものまでいくつかの式が用意されており，その中から好みのものを選ぶことができる．

図 13.8: 一般曲線座標格子

19 自明の関係 x = x( (x y) (x y)) を x で 2 階微分した式と y で 2 階微分した式
x2 x +2x x x + x2 x + xx x + xx x = 0
 2 x +2y y x +  2 x + yy x + yy x = 0yy(13.66) と (13.84) を用いれば次式が得られる．
11g x +2g12 x + g22 x + x P + x Q = 0
この式にヤコビアン J をかけ，Jg 11 = g~22  Jg 12 = ;~g12  Jg22 = g~11 なる関係を用いれば，式 (13.85a) を導出できる．同様に式
y = y((x y) (x y)) から始めれば式 (13.85b) を導出できる．
を加え，式

60

非圧縮性流れの解法|MAC 型解法

式 (13.85) の差分方程式は xy の格子点の値 xij または yij の連立 1 次方程式で，SOR 法で容易に解くことができる．ただしこれらの差分方程式の係数と右辺には解 xy が含まれるので，これらの係数と右辺の値を更新し SOR 法の計算を続行することになる．

X P U V すな
わち (2if +1)  (2jf +1) 個の点における座標 x y ，測度 x  y  x  y ，ヤコビアン J = x y ; x y ，測度
x = y =J y = ;x =J x = y =J y = x =J の値が求められる．また点 U または V における測度テンソル成分 g~11 = (x2+y2 )=J g~12 = (x x +y y )=J g~22 = (x2+y2 )=J が計算され，それぞれ配列 gU または gV に格納される．g~11 = g~22  g~12 = ;g~12  g~22 = g~11 . 最後に点 U の測度の微分 @ x =@  @ x =@ @ y =@  @ y =@
と点 V の @x =@  @x =@ @y =@  @y =@ が計算され，配列 dxix に格納される．
~ V~ と静圧 p の予測値を与えるものである．まず上流境界の流速が平均流速 1 のPREDCT は，体積流束 U
~0j =放物線分布 u0j = 6y0~j (1;y0~j ) y~0j = y0j =y0jf  v0j = 0 によって与えられる．このとき体積流速は U(J x u)0jU  V~0j = 0 となる．反変速度が  空間の流速であることを考えれば，予測値は U~ij = U~0j  V~ij = V~0jのように置くことができる．平行平板間流路の層流では，粘性による圧力降下は p = (@ 2 u=@y 2)x から求められる．ステップの上流側では @ 2 u=@y 2 = ;4=3，下流側では @ 2 u=@y 2 = ;9=16 である．また拡
大部の (回復) 圧力上昇は p = 1=2 ; :752=2 = 7=32 となる．したがって下流境界の圧力を 0 とすれば ，上流境界の圧力の予測値は p0 = f;(4=3)x0 +(9=16)xif g=Re ; 7=32 となる．プログラムでは，p の予測値
は j 方向に一様で，上流側 (0i 35) では pij = p0 +(4=3)(x0j ; xij )=Re，下流側の (35 i i0 )では pij = (9=16)(xif j ; xij )=Re とし ，またステップ下流 (35 < i < 100) では p35j と p100j から線形補R間している．なお拡大部の圧力上昇は，平均流速の u
2=2 からではなく (1=2Y ) 0Y u(y)]2 dy から求めれば ，
p = :6 ; :6  :752 = :2625 となる．
COMPUZ は点 X の流速 u v と渦度 を計算するサブルーチンである．これらの量の計算式は式 (13.68) ま
METRIC は，測度，ヤコビアン，測度テンソルなどの値を計算するものである．始めに点たは (13.72) から次のようになる．

u = y U~ ; y V~  v = ;x U~ + x V~
(13.86)no@ (~g U~ +~g V~ )(13.87)= J1 @@ (~g21 U~ +~g22V~ ) ; @1112
なお INTP は格子点 X の値が既知のときに点 U または V の値を補間するサブルーチンで，もちろん点 U ま
たは V の値から点 P の値の補間にも使用できる．INTVP は点 V の V~ の値から点 P の V~P の値を境界条件~ )i0 = 0 を考慮して補間するものである．同様に INTVX は点 V の V~ の値から点 X の V~X の値
V~i0 = (@ V=@
~ の値から点 X の U~X の値を境界条件を考慮して補間するものである．等間隔計算点
を，INTUX は点 U の U
xi の関数値 ui = u(xi ) が与えられるときに中間点 xi+1=2 の微分値は周知のように ui+1=2 = (ui+1 ;ui )=xである．また境界 x 1=2 のところでは内点の関数値を用い u 1=2 = (;2u0 +3u1 ; u2 )=x となるが，特に境界値 u 1=2 = 0 の場合には u 1=2 = (3u0 ; u1 =3)=x となる．
@ V~ 
r u = J1 @@ U~ + @
0

0
;

;

;

0
;

COMPU1 と COMPU2 は，Navier-Stokes 方程式を 2 段階に分けて解く SMAC 法において，それぞれ前段の~U = JUU ，後段の U~ = JUU を計算するサブルーチンである．前段の COMPU1 から説明する． 形陰解法の式 (13.81a) は，2 次元の場合には次のようになる．
@  ; D oU = rhsn 
J + t U~ @@ + V~ @111n@  ; D oU = rhsn
J + t U~ @@ + V~ @
2
1
2

n





61

バックステップ流路流れのプログラム
ただし

Di は粘性項の主要部で，また右辺 rhs は次のようになる．

on ~ @~+ @ V U ; u U~ @@ x + V~ @@x ; v U~ @@ y + V~ @@y + g~11 p +~g12p +  
rhs1 = ;t @@ UU
n ~ @~
x + V~ @x ; v U~ @y + V~ @y +~g 21 p +~g 22 p ; orhs2 = ;t @@ UV+ @ V V ; u U~ @



@
@
@
@

これらの式は左辺の演算子に 1 次上流差分，右辺の対流項に Chakravarthy-Osher TVD スキームを用い差分方程式に書き換えれば次のようになる．





22
22
U~ + + g~Re
(Uij ; Ui 1j ) + U~ ; g~Re
(Ui+1j ; Uij )



11
11
+ V~ + + g~Re
(Uij ; Uij 1 ) + V~ ; g~Re
(Uij+1 ; Uij ) = rhs1



22
22
JV + t U~ + + g~Re
(V ij ; V i 1j ) + U~ ; g~Re
(V i+1j ; V ij )



11
11
+ V~ + + g~Re
(V ij ; V ij 1 ) + V~ ; g~Re
(V ij+1 ; V ij ) = rhs2
JU + t





;





;





;





;





;





;






;





;

(13.88a)(13.88b)dd~ ijP ; UU
~ i 1jP + Vd
~ U ij+1X ; Vd
~ U ijX
rhs1 = ;t UU; uijU U~ @@ x + V~ @@x ijU ; vijU U~ @@ y + V~ @@y ijU11+ (~g22 )ijU (pij ; pi 1j ) ; (~g12 )ijU 2 (pij+1U ; pij 1U ) + Re ( ij+1 ; ij )
;

;

;


d
d
~ i+1jX ; UV
~ ijX + Vd
~ V ijP ; Vd
~ V ij 1P
rhs2 = ;t UV~ @x ~ @x~ @y;(13.88c)y; uijV U @ + V @ ijV ; vijV U @ + V~ @@ ijV1(; (~g12 )ijV 12 (pi+1jV ; pi 1jV ) ; (~g11 )ijV (pij ; pij 1 ) ; Re;)
(13.88d)
i+1j ij;;

c の付いている量は数値流束で，このサブルーチンには中心

なおここでは添字 ij n を一部省略している．

差分と Chakravarthy-Osher TVD スキームの２つが組み込まれている．後者の数値流束はn Ui+1=2j 1

d~ ijP = UU
~ ijP + jU~ijP j ;UU+ 6 minmod(Ui 1+1=2j  4Ui+1=2j )o+ 31 minmod(Ui+1=2j  4Ui 1+1=2j )(~Ui 1+1=2j = Ui 1=2j (U~ijP  0)
Ui+3=2j (UijP < 0)2



Ui+1=2j = Ui+1j ; Uij (13.89);のようになる．
境界条件は次のように考慮される．流路の入口と出口に逆流はないものとする．境界とその近傍点の数値

62

非圧縮性流れの解法|MAC 型解法

流束は次のようになる．

d~ 0jP = U~0jP U0j + U1j UU
2~ U i0X = Vd
~ U ijf X = 0
Vd
d
~
UV
= U~ V0j + V1j 
1jX(13.90a)Ui2 ~ U i1X = V~i1X 1 3Ui0 +2Ui1 ;
Vd
4
5 

(V~i1X > 0)d
~ i0 jX = U~i0 jX 1 (;Vi0 2j +6Vi0 1j +3Vi0 )
UV
8

1jX

2
d
~V V i0P = V~i0P 1 (12Vi1 ; Vi2 )32;(13.90b)(13.90c);(V~i0P > 0)

(13.90d)

なお最後の式は壁面上の反射の条件 @V=@ = 0 を用いて導かれたものである20~  とその隣接点の V~  の式の右辺の対流項は，非保存形の式を 1 次上流差分で近似した安下流境界点の U

定かつ境界から反射のない次式によって求められる．

~ @u + ( l )y U~ @v + V~ @v
Cl = J r l  (u ru) = ( l )x U~ @u
+
V
@
@
@
@

(l = 1 2)(13.91a)l = 1 の式の (@u=@ )if j と l = 2 の式の (@u=@ )i0 j は次のように上流差分で近似され る．(@v=@ )if j ，(@v=@ )i0 j についても同様である．
@u  = 1 (u ; u@u  = 1 (u ; u + u(13.91b)
i
i
i
0j
fj
f j +1 ; ui0 j +1 )
@ if j 2@ i0 j 2 if j i0 1j );

付加項，圧力項，粘性項の計算については説明を省略する．この計算には METRIC で計算し た dxix =@ x =@  : : : と gU,gV = g~11     が用いられる．

SMAC  形陰解法の式 (13.88) は 13.5 節に述べた改良型近似因子法を用いて効率的に解かれる．13.5 節の因子化する前の差分方程式

c0 u00 ; c+1u 10 + c1 u10 ; c+2u0 1 + c2 u01 = rhs
;

;

;

;

;
+ ;
と上記の式 (13.88a) を比較すれば ，この式の係数 c+
1  c1  c2  c2  c0 は式 (13.88a) の係数とc+1 = ;t (U~ + +~g22=Re) c1 = t (U~ ; g~22=Re)
c+2 = ;t (V~ + +~g11=Re) c2 = t (V~ ; g~11=Re) c0 = J ; c+1 ; c1 ; c+2 ; c2
;

;

;

;

;

;

のように対応していることが分かる．これらの係数は配列 c に入れられ，このようにして得られた連立 1次方程式はサブルーチン GAUSSZ で解かれる．GAUSSZ については 13.5 節の説明参照．

~
後段の COMPU2 は式 (13.81b) から U

n+1を計算する短いサブルーチンである．式 (13.81b) の差分方程式

は次のようになる．

on
U~ijn+1 = U~ij ; t (~g22 )ijU (ij ; i 1j ) ; (~g12 )ijU ij+1U ;2 ij 1Uno
V~ n+1 = V~ ; t ; (~g ) i+1jV ; i 1jV + (~g ) ( ;  )
;



ij

;

12 ijVij

2

11 ijV;~
なお非定常流れの場合にも，式 (13.83a) から U

(m)



ij

ij 1
;

~
，式 (13.83b) から U(m)

(13.92a)(13.92b)を全く同じに要領で計算することができる．

y=020
壁面上で流速成分
はそれぞれ近いところでこのような取り扱いは妥当と言える．

u vu y v y2 のようにふるまうので，曲線座標格子が壁面上で直交ないしはそれに63

バックステップ流路流れのプログラム
COEF, COMPP は

 の差分方程式の係数を設定し，その連立 1 次方程式を SOR 法または Chebyshev SLOR法で解いて  の値を求め，静圧 p の値を更新するサブルーチンである．解くべき  の微分方程式 (13.81c)
または (13.83c) は@ t g~ @ ; t g~ @  + @ ; t g~ @ + t g~ @  = @ U~ + @ V~
22
12
12
11
@
@
@ @
@
@
@
@




(13.93)のようなもので，その差分方程式は形式的に次のように表すことができる．
1
1
X
X

i = 1j = 1
0;0

cij i j i+i j+j = rhsij
0 00(13.94)0;

より具体的に書けば次のようになる．

(g1 )i+1jU (i+1j ; ij ) + (g2 )i+1jU (i+1j+1 + ij+1 ; i+1j 1 ; ij 1 ); (g1 )ijU (ij ; i 1j ) ; (g2 )ijU (ij+1 + i 1j+1 ; ij 1 ; i 1j 1 )+ (g4 )ij+1V (ij+1 ; ij ) + (g3 )ij+1V (i+1j+1 + i+1j ; i 1j+1 ; i 1j ); (g4 )ijV (ij ; ij 1 ) ; (g3 )ijV (i+1j + i+1j 1 ; i 1j ; i 1j 1 )= U~i+1j ; U~ij + V~ij+1 ; V~ij
;

;

;

;

;

;

;

;

;





;



;

;

;

;



ただし (g1 )ijU = (t g~22 )ijU  (g2 )ijU =;(t g~12 )ijU =4 (g3 )ijV = ;(t g~12 )ijV =4 (g4 )ijV = (t g~11 )ijV

で，これらの値は配列 gw に格納される．なお局所時間ステップ法では t は局所時間を取る．

 の係数 cij i j の値は，図 13.9 に示すように４段階に分けて配列 co に入力される．この図には２つの場合，領域内部の格子セル，上流境界に隣接する格子セルの場合が取り上げられている． の 1 階微分は，一般には中心差分 ( )ij = (i+1j ;i 1j )=2 が用いられるが，上流境界では片側差分 ( )0j = (;30j +41j;2j )=2が用いられる．係数 c0j 4j の値は空いている co(0,j,1,jp) に入れられる．上方と下方の流路壁に隣接す
0 0

;

0

る格子セルの場合にも，上流境界に隣接する格子セルの場合と同様の取り扱いが必要である．なお下流境
界には pif j = 0 を与えているので，下流境界に隣接する格子セルは領域内部のものと同じに取扱える．

 の連立 1 次方程式は，ここでは趣向をかえ SOR 法または Chebyshev SLOR 法で解くことにする．SOR
法では，奇数回の掃引は点 ij = 00 から i0 j0 まで，また偶数回の掃引は点 i0 ; 1 j0 ; 1 から 11 まで逆向きに行われる．境界点の点緩和は，このように過緩和の効果が損なわれないように行われるべきである． の差分方程式 (13.94) の残差値はResij =

1
1
X
X

i = 1j = 1
0;0

cij i j i+i j+j ; rhsij
0 0

0

0

;

から求められ，点過緩和の修正計算は次式によって行われる．

(ijn+1) = (ijn) ; Resij =cij00
 = 1:5 は過緩和係数である．

非定常流れの場合には，各時間ステップごとに， の差分方程式を十分に満足させる必要があり，ここではベクトルコンピュータの性能を十分に引き出せるよう，Chebyshev SLOR 法を用意した．この方法では緩和計算は奇数列 (i=1,3,5, : : : ) と偶数列 (i=0,2,4, : : : ) の２段階に分け，各段階の計算では３重対角行列64

非圧縮性流れの解法|MAC 型解法
3rd step
g4
;g3

;g4

g3
g3

;

01

11

;

00

10

i

i +1

;g4

;g3

g4 ;3g3 4g3;g4 ;3g3 4g3

;g3

;g3

2nd step;g2 ;g2

g1
g2

g2

;

;

;

4th step

g3
g3

内 部の格子セル

g4

;g2

0 0

;

;

g2
g1
;g2

;g3

3rd step

j +1 0 1
j 00
j 1 i 0j =1
i=1

g2

;g1

0 0

;

;g3
1st step
11

11

;

10

10

;

1 ;1
2

1 ;1
3

;

4th step;g4 +3g3 ;4g3

g4 +3g3 ;4g3

上流境界の格子セル

g34
 g4

1st step

j +1 1 1j
10ijj 1 1=1i 1;g1

g1  g2 

g2
;g1
;g2

g2
g1
;g2

g3
g3

図 13.9: 圧力差分方程式の係数の設定の連立１次方程式の系が Gauss 消去法で解かれる．１つの列の連立１次方程式は

ci0 00 (i0n+1) + ci0 01 (i1n+1) + ci0 02 (i2n+1)= ri0 ;

2
Xn) ) ; (1 ; )c (n)(ci0 1j (in)1j + ci0 1j (i+1
i0 00 i0
j
;

0

;

0

0

0

j =0
(n+1)
cij 0 1 ij 1 + cij 00 (ijn+1) + cij 01 (ijn+1)+1
0;= rij ;

;

1
X

n)
(n)
(cij 1j (in)1j+j + cij 1j (i+1j +j ) ; (1 ; )cij 00 ij
;

0

0

0

;

j= 1
(n+1)
(n+1)
cij0 0 2 (ijn+1)
0 2 + cij0 0 1 ij0 1 + cij0 00 ij0
0;= rij0 ;0(j = 1 2 : : :  j0 )

;

;

;

X
0

n)
(n)
(cij0 1j (in)1j0 +j + cij0 1j (i+1j0 +j ) ; (1 ; )cij0 00 ij0

j= 2
0

;

;

0

;

0

0

0

;

のようになる．ただし  は加速のパラメータである．なお上流境界の連立１次方程式は境界条件が入り多少異なバックステップ流路流れのプログラム65

図 13.10: バックステップ流路の流れ (Re = 100)
るものになる．係数 cij i0 j 0 の値は一般に配列要素 co(i,j,ip,jp) に入れられるが，ci0 i0 2 は co(i,0,ip,-1)に，cij0 i0 ;2 は co(i,j0,ip,1) に入れられる．

これらの連立１次方程式はサブルーチン GAUSSP を引用して解かれる．GAUSSP では奇数列または偶数列の連立 1 次方程式が一まとめに同時に解かれる．また各連立１次方程式は，次のように係数行列が変形３重対角行列になっている．

66

非圧縮性流れの解法|MAC 型解法

図 13.11: バックステップ流路の流れ (Re = 500)

0a a a10 x 1 0 b 1
0i02 i03 i01
B
CC BB xii01 CC BB bii01 CCa
Bi
11 ai12 ai13
B
CC BB x CC BB b CC
B
ai21 ai22 ai23
B
CC BB .i2 CC = BB i.2 CC
B
.
.
.
..
..
..
B
CC BB .. CC BB .. CCB
B@
A B@xi n 1CA B@bi n 1CA
ai n 1 1 ai n 1 2 ai n 1 3 C
0

;

ain3

;

ain1

;

ain2

;

xin

;

bin

ただし添字 i は i 番目の連立１次方程式を意味する．この係数行列は３重対角要素に要素 ai01 と ain3 が加わったもので，その Gauss 消去法の計算では，通常の３重対角要素の計算に，その開始時，折返し点，終了時にこれらの要素の計算が追加される．このプログラムは，定常低レ イノルズ数の場合には SOR 法，その他の場合には Chebyshev SLOR 法で計算するように組まれている．これらの反復計算の掃引数は残差値

バックステップ流路流れのプログラム67

図 13.12: バックステップ流路の流れ (Re = 1000)
resNS と resP の大きさが同等になるように決められたものである．なお Chebyshev SLOR 法はベクトル計算機を利用するときに威力を発揮する．

次に計算結果について述べる．ここにはレ イノルズ数 100, 500, 1000 の結果を示すが ，これらはいずれ

も t = 1:0，mode = 2 と置いて計算したものである．サブルーチン COMPU1 では Chakravarthy-Osher 型の TVD スキームを用い= 1:0， = 0:5 と置き，またサブルーチン COMPP では SOR 法を用いる場合には  = 1:5， = 0:5，チェビシェフ SLOR 法の場合にはレ イノルズ数 100 では  =  = 1:0，レ イノル

ズ数 500 と 1000 では  = 1:25  = 0:75 としている．図 13.10，13.11，13.12 は SOR 法の場合の結果である．これらの図には，上から圧力分布と速度ベクトルプロフィル，渦度分布と流線が示されている．計算は x =;10  70:91 の範囲で行われたが，図には全領域ではなく x = ;6  18 と x = 18  42 の範囲が示されている．圧力分布は例えば Re = 100 の場合には，一つの色の範囲は 0.015 で，p = 0:21 の色は

p = 0:2025  0:2175 の領域に塗られている．流線は下壁面で流れ関数  = 0 上壁面で 3 とし ， = 0:25

68

非圧縮性流れの解法|MAC 型解法

おきのものが描かれている．また特に剥離域内では  = 0:01 おきに描かれている．収束判定は前章同様

resNS < 0:0001 かつ resP < 0:0001 で行っている．Re = 1000 の場合の解の収束状況を下表に示す．
表 13.6: バックステップ流路流れの解の収束状況 (Re = 1000)
10;5
SOR

Chebyshev SLOR

n

resNS

resP

flow

resNS

resP

flow

100
1000
2000
2340
2460

736
41
12
9

436
50
16
10

25693
5095
3435
2608

481
157
31

189
79
25

29678
;6564
4669

9

10

4402

収束解を得るまでの反復数 n は， 形陰解法の適用によって大幅に減少しているが，時間間隔 t やプロ

グラム中の係数  などにもより，それらの値も適切に設定しないと効率よく収束解を得ることはできない．
次表は SOR 法を用いた場合の時間間隔と収束に要する反復数の関係を示しものである．
表 13.7: 時間間隔 t と収束反復数 n

Re
100

500

1000

t

n

xrap

pmax

pmin

0.1
0.25
0.63
1.0
2.5
4.0
0.25
0.4
0.63
1.0
1.6
2.5
4.0
0.63
1.0
1.6
2.5
4.0

2580
1180
900
1366
3562
*
**
3100
1580
1744
2713
4156
*
**
2348
2917
3951
*

7.2
7.2
7.1
7.2
7.2
7.2

0.245
0.240
0.242
0.242
0.242
0.243

0.064
0.058
0.064
0.062
0.062
0.062

21.8
21.7
23.5
23.5
23.3
23.7

;0:017;0:018;0:016;0:016;0:016;0:015;0:159;0:159;0:155;0:155;0:155;0:151

30.7
38.3
38.8
41.1

;0:038;0:033;0:033;0:034;0:130;0:129;0:130;0:125

表中の*印は n = 5000 回まで計算した未収束の結果，また**印は残差値の状況から収束解が得られないと判断されたことを示す．この表にはまた再付着距離，図示領域の最大最小圧力値も示してある．これらは本
来時間間隔 t には無関係であるべきものであるが，表に示すように必ずしも無関係とはいえない．次表はこれらの値を更に，t = 1:0 に対し ，本章の Chebyshev SLOR 法と前章の長方形格子 (RG) の計算結果と比較したものである．
このことについて多くの読者は収束判定条件があまいのではと思われるかもしれない．しかし 結論を先に言えば判定条件は問題ないが ，流路の角近傍の精度が不十分ということである．このプログラムでは比較的粗い格子に高次スキームを用い精度を確保しようとしているが，高次スキームも角近傍の極端に歪んだ格子に対しては精度を発揮できない．この曲線座標格子では，上記の判定条件で計算が収束した後に更に計算を継続すると，ある時点で残差が突然増大しその後再付着距離のかなり長い解に収束するかまたは~360 の符号が正から負へ切替わることによる．通常発散するという不都合が生じる．その原因はおそらく U

69

反変速度を用いた拡散項の簡潔な表示

表 13.8: バックステップ流路流れの再付着距離，最大最小圧力の比較
RG

SOR

Chebyshev SLOR

Re

xrap

xrap

pmax

pmin

xrap

pmax

pmin

100
500
1000

7.3
20.8

7.2
23.5
30.7

0.242
;0:016;0:0380.062;0:155;0:130

7.1
22.5
31.4

0.242
;0:013;0:0330.064;0:143;0:124点では 13.4 節に詳しく述べたようにこの符号の切替わりにより対流項の値が急変することはないが，角の~360 の対流項を 2 次の TVD 補間式を微直後のこの点ではするということである．上記のプログラムでは U分した式を用いて計算するようにし ，とりあえずこの問題に対処している．しかし正攻法はより細かい格子を用いることであろう．

13.9

反変速度を用いた拡散項の簡潔な表示

初めに層流の計算や乱流の直接シミュレーションに用いられる粘性項

らJ

r2u を考える．この粘性項に左かr l  を演算したものは，式 (13.70){(13.74) を用いれば次のようになる21 ．


Dl = J r l  r2u = ; J r l r = ; lij @@ @@xx  = ; lij @@ (~gjk Z~k )iji(13.95)~k は式 (13.72) から求められる．すなわち
これは６項からなる簡潔な式である．渦度 Z

Z~k = kmn @ @ (~gni U~i )m

なお反変渦度 Zk は立方体セルの面に垂直な成分k 軸に平行な稜の中心点に定義される．また２次元の場合には渦度は紙のみとなり，粘性項は次のようになる．

@  D = @D1 = ; @2
n @ ~ ~@ @ ~ ~ o1= J @ (~g21 U +~g22V ) ; @ (~g11 U +~g12V )(13.96)

次に乱流の計算に用いられるより一般的な粘性拡散項を考える．この粘性項は粘性応力テンソル  を用

い一般に次のように表わされる22 ．

@ l @mn = @ J @ i @ l   ; J @ i  @ @ l
Dl = J r l  (r ) = J @x
@ i @xm @xn mn@xm mn @ i @xn
n @xm= @@ (JT il ) + ilj JT iji
ただし mn は  の成分，T ij はその反変成分で次のように表される．
@i @jT ij = @x
mn @x
m
n

(13.97)

21 第 2 式から第 3 式へは r(ra ) = r(ra );r2 a なる関係，第 3 式から第 4 式へは式 (13.65) と (13.71)，また第 4 式から

第5 式へは式 (13.73) と (13.74) が用いられた．

22 式 (13.64)，(13.65)，(13.78) を用いる．

70

非圧縮性流れの解法|MAC 型解法

乱流を渦粘性近似を置いて計算する場合には，応力テンソル成分は次式で与えられる．

@un + @um  ; 2  k(13.98)mn = e @x
3 mn
m @xn
ただし e = + t ， t は渦粘性係数，k は乱れの運動エネルギーである．このとき応力テンソルの反変成
分 JT ij は次のようになる．

"@ i n @ u @ j ;u @2 j on @x @x@xm @xm n @xn
n mn@o 2 ij #
2
@
@
@
j
i
i
+ e @x @x um @x ; um @x @x ; 3 g k
n
n
m
m n

@
@U
@
@
@ j @Ui ; U @ @ i  ; 2 g~ij k= e J @x i @x j ; Uk @ @x j + e J @x
k @ @x
3
m
m
k m
n @xn
k n
2
 e (Aij + Aji ) ; 3 g~ij k(i j = 1 2 3)

JT ij = Je

(13.99)ここに@ i @Uj ; U @ @ j 
Aij = J @x
k @ @xn @xn
k n @i j@@U@b= J @xi @xm @ j + Uk @@xm mj k  @x
n
n m
n
n n

= g~im @@Uj + mj k Uk
 g~im ajmm

拡散項は，渦粘性近似のもとでは，式 (13.97) と式 (13.99) から次のようになる23 ．

nonoDl = @@ e (Ail + Ali ) ; 32 g~il k + ilj e (Aij + Aji ) ; 23 g~ij ki(l = 1 2 3)= e @@ (Ail ; Ali ) + @@ e (Ail + Ali ) ; 23 g~li @@kiii(13.100)

ただし@ i bjAij = g~im ajm = @xn n
ajm = @@Uj + mj k Uk m@ k U~ ; @ j U~ 
bjn = @@ @x
j @x kknn(i j = 1 2 3)(13.101)~1 と同じ立方体セルの面中心に，A12  A21 は Z3 と同じ陵の中心に定義さである24 ．なお A11  a1m  b1n は Uれる．

23 次のように計算する．

 @ @l @@i l l
i j Aij = @xn i j bjn = ; @j @xn bjn = ; @j Alj @i @ @l@li j g~ij = ;J @xn @i @xn = ; @i g~il

24 @bjn =@j = 0，したがって bjn は j = const: 線に沿って一定になる量である．