第14 章 非圧縮性流れの解法|MAC 型解法(続き)14.1 直角座標系から一般曲線座標系への変換
直角座標 (cartesian coordinates) を x ，一般曲線座標 (general curvilinear coordinates) を  で表す1 ．これらの座標間の変換の測度 (metric) x  : : : と x  : : : の間には次の関係が成立する2 ．

2
x
66
4y

z

32

3

2

3

x x 7 6x y z 7 6 1 0 0 7y y 75 64x y z 75 = 64 0 1 0 75z z x y z
0 0 1(14.1)

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

x = Jy  x = ;Jy  y = ;Jx  y = Jx J = x y ; x y

x に関する微分を  に関する微分に書き換えれば次のようになる3 ．
@ = @j @  r = (r ) @j @
@xi @xi @jj
具体的には，例えば y@ = 1  z
@x J 




y
y 

z
z  @=@ @=@ @=@ (14.3)(14.4)@ = 1 y @ ; y @@x J  @ @

２次元では

また一般に

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

非圧縮性流れの解法|MAC 型解法であるから次の関係が導かれる．

@ = J @j @ = @ J @j
J @x@xi @j @j @xii@J r = @ J (rj )j@2
J r = @ Jgjk @@j(14.5a)(14.5b)(14.5c)k

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

@i @jgij = ri rj = @x@xk(14.6)

k

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

@j @ = U @ = U  r~
u r = ui @x@ = ui @x
j @
@
i

i

j

(14.7)j

ただし Uj は反変速度 (contravariant velocity) の j 方向成分である．反変速度 U は  空間の流速で 4 ，物理空間の流速 u との間には次の関係がある．

@j u Uj = @xiiui = @x@ Uji(14.8)

j

なおこのような関係を知らずに，基礎方程式の xi の微分を，式 (14.4) によって j の微分に書き換えるとたいへん長い式になり収拾がつかなくなる．
流れの渦度  =ru の xi 方向成分はk
i = ijk @u@x(14.9)j

ただし ( ijk ) は Eddington の で

8>>< 1 ((ijk) = (123) (231) または (312))
ijk = > ;1 ((ijk ) = (321) (213) または (132))>: 0 (その他の場合)

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

; @u2=@x3 である．任意のベクトル a に対し次の関係が成立する．
@xi @ai = 1@ @xx arl ra = rl rm  @@aa = J1 lmn @
lmn@m @nm
n @m J

(14.10)

なおこの式の導出には次式が用いられた．

y z x = J   z z @ 2x

lmn @ @ = 0
m n(14.11)式 (14.10) で，特に a = u ならば ．次の反変渦度と反変速度の関係が得られる．

Zl = J1 lmn @@ (gnj Uj )(14.12)

m

x

x +u

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

x



x +u

 +u 

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

3

ただし (gij ) は測度テンソルで，(gij );1 = (g ij )，その成分は次式で定義される．

@xx  @xx = @xk @xk
gij = @@ @ @iji

(14.13)j

また Zl は反変渦度の l 方向成分で，反変渦度 Z と物理空間の渦度  の間には次の関係がある．

@j  Zj = @xiii = @x@ Zji(14.14)j14.2 曲線座標格子の SMAC 形陰解法
物理空間の曲線座標格子を写像空間の立方体格子に写像する．立方体格子の稜の長さ   は通常 1 に取られるが，以下の説明では当面 省略せずに残しておくことにする．直角座標格子の MAC 型解法

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

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

i = @ J @j u = @ U~ = 0
J @u@x @ @x i @ ji

ji

(14.15)j

図 14.1 は  空間内の１つの立方体格子セルを示したものである．式 (14.15) は，この立方体格子セルにわたって積分すれば(U~100 ; U~000)   + (V~010 ; V~000 )   + (W~ 001 ; W~ 000 )   = 0(14.16)となる．式 (14.16) の中で， u x x @ u
U~000   = J @x
  =    v y y ii 000Uw z z  u000U  (x 010 ;x000 )  (x 001 ;x000 )000U;x 000 )(x 001 ;x 000 ) は x 空間内の対応する六面体格子セルの， = 0 面の面積ベクトルであ~るから，U000   はこの面を通る流量である．したがって式 (14.16) はこの六面体セルの６つの面から出入りする流量の和が 0 になるという意味の式である．写像空間内の流れに垂直に取られた単位面積当たり~ V~  W~ はこの体積流束の    方向成分で，の物理空間の体積流量は体積流束 (volume ux) と呼ばれ，U
また  =  =  = 1 ならば対応する六面体セルの面を通る流量になる．なお上式の導出には式 (14.1)

ここに (x 010から得られる次の関係が用いられた．

y y 1x = J  z z 4

非圧縮性流れの解法|MAC 型解法~l は流速 u の左から J
次に体積流束の運動方程式について述べる．体積流束成分 Url  を演算したもので~l の運動方程式は運動方程式 duu=dt  @uu=@t+r  uu = ;rp+ r u +g の左から J rl  を演算あるから，U
2

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

Cl = J rl  (ruu) = rl  @@ U~iu = @@ U~i Ul ; u  U~i @@ rl
i
i
i


@
l
= @ U~i Ul + i j 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(14.17)(14.18)なお式 (14.17) の導出には次の関係が用いられた6 ．

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

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

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

@p ; D n
U~l = U~ln ; t Cl +~gli @li@
U~ln+1 = U~l ; t g~li @i@ g~li @ = 1 @ U~   pn+1 = pn +@@t @ llil(l = 1 2 3)(14.20a)(14.20b)(14.20c)x001
HH 6HHHHH
W~ 
001
U~H
000H
H


:

HHj V~010
6

:

:

U~100

V~000 H
HHj
HH
x

010 H
j H
6
H
x000 H
HHHW~ 000Hx100
図 14.1:  空間内の立方体格子セル~i Ul = l (U~iu) の微分が用いられた．
2 式から第 3 式へは式 (14.5) と (14.8)，第 3 式へは U

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

r

r

r

5

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

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

noJ + t U~i @@ ; Dl Ul = rhsnl (l = 1 2 3)i@p ; D  U~  = U~ n + J U rhsl = ; t Cl +~gli @
l
l
l
l
i
@
U~ln+1 = U~l ; t g~li @i@ g~li @ = 1 @ U~   pn+1 = pn +@@t @ lli(14.21a)(14.21b)(14.21c)l~l = U~ln であるから，式 (14.21a) の右辺 rhs は残差，定常流れの計算では，解が収束すれば U

Ul は修正値を表している．解が収束すればこの残差はゼロに近づき，修正値もゼロに近づく．したがってこの式の右辺は，解の精度に直接関係するので精度良く離散化しなければならないが，左辺の演算子は大胆に近似することができ，1 次上流差分が用いられ，拡散項は主要部のみ考慮され，また近似因子法が適用される．

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

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

n 


o
D1 = e @@ g~22 g~33 @@ J  + @@ g~33 @@ g11 + @@ g~22 @@ g11(14.22)また２次元の場合にはn o2
D1 = 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) (14.23a)
i

(
m
)
(
m
;
1)
U~l = U~l
+ J Ul(m)(m)
U~l(m) = U~l(m) ; 21 t g~li @ @
(14.23b)i@ g~li @ (m) = 2 @ U~ (m)  p(m) = p(m;1) + (m)(14.23c)l
非定常流れの計算では，

@l

@i

t @l

~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 時間ステップ当たりの反復回数は少しで済むことになる．

6

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

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

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


START 
GRID

METRIC
COEF

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

m = m +1e

COMPU1COMPPCOMPU2COMPUZ

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

f

CG


STOP 図 14.2: SMAC  形陰解法のフローチャート

7

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

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

~ ，COMPP ででU~ ，また COMPUZ で  の計算が行われる．n20 回ごとに NS 差分方程式と p，COMPU2 で Uの最大残差値と圧力差分方程式の最大残差値を打出し ，また非定常流れでは CG に必要なデータをファイルに書込む．この反復計算は解が発散したときには直ちに中止し ，解が収束または規定の反復数に達した
時に次のステップまたは計算を終了する．終了直前には計算の最終結果をファイルに書込み CG で流れの可視化を行う．以下にバックステップ流路を通る定常流れのプログラムを示す．
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! *****************************************************************************************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
CALL COMPUZ(UJ,u,z,locf,fmax,CFLmax)!undteady: initial data!velocities vorticity & divergence

OPEN(20,FILE='OUTPUT.dat')
WRITE(20,60) IF(unsteady)WRITE(20,61)n=0 100 n=n+1!n: stepm=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)
CALL COMPUZ(UJ,u,z,locf,fmax,CFLmax)! Decide convergence and output computational resultsIF(MOD(n,20)==0 .AND. (steady.OR.m==1))THEN
i=if resNS=0.

!compute JU!velocities vorticity & divergence8

非圧縮性流れの解法|MAC 型解法
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!convergeIF(steady)

GOTO 110IF(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= 0

CASE(2) FORALL(i=0:if,j=0:jf)f(i,j)=x(2,i,j)  z1="y" z2="" mz= 0

CASE(3) FORALL(i=0:if,j=0:jf)f(i,j)=u(1,i,j)  z1="u" z2="" mz= 0

CASE(4) FORALL(i=0:if,j=0:jf)f(i,j)=u(2,i,j)  z1="v" z2="" mz= 0

CASE(5) FORALL(i=0:if,j=0:jf)f(i,j)=p(i,j)

 z1="p" z2="*10 " mz= 1

CASE(6) FORALL(i=0:if,j=0:jf)f(i,j)=z(i,j)

 z1="zeta" z2="/10 " mz=-1

ENDSELECTCALL WRTF(f,z1,z2,mz,n)!mz:scale factor!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

9

バックステップ流路流れのプログラム
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

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_xi, y_xi
!x_eta, y_eta

FORALL(j=0:jf)xxi(l, 0,j)=-(3.*x(l, 0,j)-4.*x(l, 1,j)+x(l,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

FORALL(i=0:if)xxi(l+2,i, 0)=-(3.*x(l,i, 0)-4.*x(l,i, 1)+x(l,i,2,j))/2./dxi2))/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

10

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

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/67/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/69/6)!

Q(i,j)=coef*(18.*ty-10.)!te(ty)=3ty^3-5ty^2+3ty

te'=(18/6 1.5/6 12/6)

4/6
3/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))
i=if DO 13 j=1,j0!outlet boundw=(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))
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)))!|dx|+|dy|

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
! ********** 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), &11

バックステップ流路流れのプログラム
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
xix(1,i,j)= xxi(4,i,j)/aJ(i,j) xix(2,i,j)=-xxi(3,i,j)/aJ(i,j)!xi_x, xi_y at X

xix(3,i,j)=-xxi(2,i,j)/aJ(i,j) xix(4,i,j)= xxi(1,i,j)/aJ(i,j)!eta_x, eta_y at X

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)!(xi_x)_xi

at U

dxix(i,j,1,3)=xix(2,ii+1,jj)-xix(2,ii-1,jj)!(xi_y)_xi

at U

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)!(xi_x)_eta at U

12非圧縮性流れの解法|MAC 型解法
dxix(i,j,1,4)=xix(2,ii,jj+1)-xix(2,ii,jj-1)!(xi_y)_eta at U

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)!(eta_x)_xi

at V

dxix(i,j,2,3)=xix(4,ii+1,jj)-xix(4,ii-1,jj)!(eta_y)_xi

at V

dxix(i,j,2,2)=xix(3,ii,jj+1)-xix(3,ii,jj-1)!(eat_x)_eta at V

dxix(i,j,2,4)=xix(4,ii,jj+1)-xix(4,ii,jj-1)!(eta_y)_eta at V

ENDDO ENDDO
ENDSUBROUTINE METRIC! ********** 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
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)

13

バックステップ流路流れのプログラム
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

FORALL(i=0:if)w1(i)=UJ(1,i,j) CALL INTP(w1,w11,if)!w11=JU_P=U

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.)*duELSE 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

14

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

i=0  hf(i)=0.
i=1  hf(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

FORALL(j=0:jf)w2(j)=UJ(2,i,j) CALL INTVP(w2,w21,jf)!w21=JV_P=V

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)!residuals of NS eqn= Uq(1,i,j)-ad+pr+(z(i,j+1)-z(i,j))/Re

rhs(1,i,j) = -dt*Uq(1,i,j)
ENDDO ENDDO!rhs of linear eqns

15

バックステップ流路流れのプログラム
DO j=1,j0
FORALL(i=0:if)w1(i)=u(1,i,j)

CALL INTP(w1,w11,if)!w11=u_V

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)!residuals of NS eqn= Uq(2,i,j)-ad+pr-(z(i+1,j)-z(i,j))/Re

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

16

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

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
co(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)!1st step

17

バックステップ流路流れのプログラム
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)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)!bottomco(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)!top

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!3rd stepco(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)!inletco(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!4th stepco(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)!inlet18

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

ENDDO
ENDSUBROUTINE COEF
! ********** 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
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)-cop19

バックステップ流路流れのプログラム
r(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
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)! ***** Compute odd j-rows
DO i=0,i0
DO jj=1,jfo j=2*jj-1!even row phi

20

非圧縮性流れの解法|MAC 型解法
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)IF(resP>1.E-5 .AND. n<nf)!odd row phi
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), &21

バックステップ流路流れのプログラム
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=0  z(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

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
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)RETURN!max|CFL|

22

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

IF(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)
IF(ABS(fa)<.00001)!fa=f(alpha) by cubic interpolation
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
b = b-fb/dfb

!dfb=f'(beta)
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,ifバックステップ流路流れのプログラムb(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=0

 u1(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)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)

23

24

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

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=0

 u1(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(","DATA area/" -3P"," -2P"," -1P","","","F8.4","))
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),
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

&

25

バックステップ流路流れのプログラム
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!greencolors(6) =#60FFA0 colors(7) =#50FFD0 colors(8) =#40FFFF!yellowcolors(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

26

非圧縮性流れの解法|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)
CALL LINE(-6.0,0.0,0.0,0.0)!gray 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

27

28

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

numfonts=INITIALIZEFONTS()
fontnum=SETFONT("t'Arial''h16w6'")
xt=x yt=y+0.33 CALL MOVETO_W(xt,yt,wxy) CALL OUTGTEXT(A)ENDSUBROUTINE CHARAC! ********** 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)

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
i1=i!for all elements

j1=j

x1=x(i1,j1) y1=y(i1,j1) f1=f(i1,j1)

i2=i+1 j2=j

x2=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+1m(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)!for all quadrilateral elementsm(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.

29

30

非圧縮性流れの解法|MAC 型解法! 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!

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
xs=(a1-a)/(a1-a4)*x4+(a-a4)/(a1-a4)*x1!ntype=2!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
yt=(a1-a)/(a1-a2)*y2+(a-a2)/(a1-a2)*y1
ELSEIF(k<k3)THEN!k1<=k<k2バックステップ流路流れのプログラム
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!k==k1<k2

IF(k==k2)THEN
poly(4)%wx=x2 poly(4)%wy=y2 nvertx=4!k==k1==k2<k3

ENDIF
IF(k==k4)THEN!k==k1==k2==k3==k4

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!k==k1==k2==k3<k4

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!k1<k==k431

32

非圧縮性流れの解法|MAC 型解法
poly(3)%wx=x4

poly(3)%wy=y4 nvertx=3

IF(k==k3)THEN

nvertx=4!k3<k==k4!k2<k==k3==k4

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!k1<k==k2==k3==k4

poly(5)%wx=x2 poly(5)%wy=y2 nvertx=5
ENDIF
ELSE

!k1<k<k4poly(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!k/=k2,k3

IF(k==k2.AND.k<k3)THEN
poly(5)%wx=x2 poly(5)%wy=y2 nvertx=5!k1<k==k2<k3

ENDIF
IF(k==k3)THEN

nvertx=5!k2<k==k3<k4

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
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!k1<k==k2==k3<k433

バックステップ流路流れのプログラム! ********** 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)
ENDSUBROUTINE LINE

status = LINETO_W(xt,yt)

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

140 0  25 に写像される．格子間隔  =  = 1．上流境界 (x = ;10) の格子点番号
は i = 0，下流境界 (x = 70) は if = 140，ステップ (x = 0) の角は i1 = 35，隅は i2 = 40 である．また下方壁面 (y = 0 (x 0) x = 0 (;1 < y < 0) y = ;1 (x  0)) の格子点番号は j = 0，上方壁面 (y = 3) は
jf = 25 である．
曲線座標格子は基本的に次の Poisson 方程式の境界値問題を解くことによって形成される．
領域 0xx + yy = P

xx + yy = Q(14.24)

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

ただし L は微分演算子で
2
2@2
L = g~22 @@ 2 ; 2~g12 @ @@ +~g11 @
2g~ij = gij =J ，ヤコビアン J と測度テンソル gij は前節で定義されたものである7 ．

式 (14.25) の境界値問題の境界条件は次のように与えられる．下方境界には境界点の座標が等比級数に;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

7 自明の関係 x = x( (x y) (x y)) を x で 2 階微分した式と y で 2 階微分した式

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

34

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

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

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

i 10
;

1
;x ; x ;x1

xi+10

i0

i0



i 10
;

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

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

00

ij

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

0

0

おプログラムには，格子間隔の差の小さいものから大きいものまでいくつかの式が用意されており，その中から好みのものを選ぶことができる．
式 (14.25) の差分方程式は xy の格子点の値 xij または yij の連立 1 次方程式で，SOR 法で容易に解くことができる．ただしこれらの差分方程式の係数と右辺には解 xy が含まれるので，これらの係数と右辺の値を更新し SOR 法の計算を続行することになる．

図 14.3: 一般曲線座標格子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
METRIC は，測度，ヤコビアン，測度テンソルなどの値を計算するものである．始めに点35

バックステップ流路流れのプログラムのように置くことができる．平行平板間流路の層流では，粘性による圧力降下は p =(@ 2 u=@y2)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 の予測値
i 35) では pij = p0 +(4=3)(x0j ; xij )=Re，下流側の (35 i i0 )は j 方向に一様で，上流側 (0では pij = (9=16)(xif j ; xij )=Re とし ，またステップ下流 (35 < i < 100) では p35j と p100j から線形補R
2=2 からではなく (1=2Y ) 0Y u(y)]2 dy から求めれば ，間している．なお拡大部の圧力上昇は，平均流速の u
p = :6 ; :6  :752 = :2625 となる．
COMPUZ は点 X の流速 u v と渦度 を計算するサブルーチンである．これらの量の計算式は式 (14.8) または (14.12) から次のようになる．

u = y U~ ; y V~  v = ;x U~ + x V~
(14.26)no@ (~g U~ +~g V~ )(14.27)= 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 となる．
0

0
;

;

;

@ V~
r u = J1 @@ U~ + @

0;

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

o
@ ; D U = rhsn
J + t U~ @@ + V~ @
2
1
2




ただし

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

n ~ @~
orhs1 = ;t @@ UU+ @ V U ; u U~ @@ x + V~ @@x ; v U~ @@ y + V~ @@y + g~11 p +~g12p +  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 スキームを用い差36

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

分方程式に書き換えれば次のようになる．





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


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


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


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




;





;





;





;





;





;






;





;


d
d
~ ijP ; UU
~ i 1jP + Vd
~ U ij+1X ; Vd
~ U ijX
rhs1 = ;t UU~@ y~@ x ~@ x;(14.28a)(14.28b); uijU U @ + V @ ijU ; vijU U @ + V~ @@y ijU1(;)+ (~g22 )ijU (pij ; pi 1j ) ; (~g12 )ijU 21 (pij+1U ; pij 1U ) + Reij +1 ij;;(14.28c)dd~ i+1jX ; UV
~ ijX + Vd
~ V ijP ; Vd
~ V ij 1P
rhs2 = ;t UV@y + V~ @y
x + V~ @x ; v~; uijV U~ @U@@ ijV ijV @@ ijV11; (~g12 )ijV 2 (pi+1jV ; pi 1jV ) ; (~g11 )ijV (pij ; pij 1 ) ; Re ( i+1j ; ij ) (14.28d);;;

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=2j (U~ijP  0)
Ui 1+1=2j =Ui+3=2j (U~ijP < 0)2



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

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

(14.30a)(V~i1X > 0)d
~ i0 jX = U~i0 jX 1 (;Vi0 2j +6Vi0 1j +3Vi0 )
UV
2
8
1
d
V~ V i0P = V~i0P 32 (12Vi1 ; Vi2 )(V~i0P > 0)
1jX

1jX

;

;

(14.30b)(14.30c)(14.30d)

なお最後の式は壁面上の反射の条件 @V=@ = 0 を用いて導かれたものである88 y = 0 壁面上で流速成分 u v はそれぞれ u  y v  y2 のようにふるまうので，曲線座標格子が壁面上で直交ないしはそれに近いところでこのような取り扱いは妥当と言える．

37

バックステップ流路流れのプログラム~  とその隣接点の V~  の式の右辺の対流項は，非保存形の式を 1 次上流差分で近似した安下流境界点の U
定かつ境界から反射のない次式によって求められる．

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

(l = 1 2)(14.31a)l = 1 の式の (@u=@ )if j と l = 2 の式の (@u=@ )i0 j は次のように上流差分で近似され る．(@v=@ )if j ，(@v=@ )i0 j についても同様である．
@u  = 1 (u ; u + u@u  = 1 (u ; u(14.31b)
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  形陰解法の式 (14.28) は 13.5 節に述べた改良型近似因子法を用いて効率的に解かれる．13.5 節の因子化する前の差分方程式

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

;

;

;

;
+ ;
と上記の式 (14.28a) を比較すれば ，この式の係数 c+
1  c1  c2  c2  c0 は式 (14.28a) の係数と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 は式 (14.21b) から U

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

は次のようになる．

noU~ijn+1 = U~ij ; t (~g22 )ijU (ij ; i 1j ) ; (~g12 )ijU ij+1U ;2 ij 1Uon
V~ n+1 = V~ ; t ; (~g ) i+1jV ; i 1jV + (~g ) ( ;  );12 ijVij

ij(14.32a)

;





2

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

(m)



ij

(14.32b)

ij 1
;

~
，式 (14.23b) から U(m)を全く同じに要領で計算することができる．

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



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

i = 1j = 1
0

;

0

;

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

0

0

(14.34)38

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

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

(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 は局所時間を取る．
3rd step
g4
;g3
;g3

2nd step;g2 ;g2

g1
g2

;g4

g2

;

01

11

;

00

10

i

i +1

;g4

;g3

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

;g3

;

;

;

4th step

g3
g3

g4

3rd step

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

0 0

;

10
1 ;1
2

;g4 +3g3 ;4g3

g4 +3g3 ;4g3

g2
;g1
;g2

g2
g1
;g2

;g3

;g3
1st step11

4th step

上流境界の格子セル

g34
 g4

0 0

;

内 部の格子セル

g1  g2 

1st step

j +1 1 1j
10ijj 1 1=1
i 1

;g1

g3
g3

11

;

10

;

1 ;1
3

;

g2
;g1
;g2

g2
g1
;g2

g3
g3

図 14.4: 圧力差分方程式の係数の設定

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

;

39

バックステップ流路流れのプログラムが用いられる．係数 c0j 4j 0 の値は空いている co(0,j,1,jp) に入れられる．上方と下方の流路壁に隣接する格子セルの場合にも，上流境界に隣接する格子セルの場合と同様の取り扱いが必要である．なお下流境
界には pif j = 0 を与えているので，下流境界に隣接する格子セルは領域内部のものと同じに取扱える．

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

Resij =

i = 1j = 1
00;

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, : : : ) の２段階に分け，各段階の計算では３重対角行列の連立１次方程式の系が 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
cij 0 1 (ijn+1)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
cij0 0 2 (ijn+1)+ cij0 0 1 (ijn+1)+ cij0 00 (ijn+1)
0 2
0 1
0
0

;

= rij0 ;0(j = 1 2 : : :  j0 )

;

;

;

0
X

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

;

のようになる．ただし  は加速のパラメータである．なお上流境界の連立１次方程式は境界条件が入り多少異なるものになる．係数 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 次方程式が一まとめに同時に解かれる．また各連立１次方程式は，次のように係数行列が変形３40

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

重対角行列になっている．

0
10 1 0 1a
0i
02 ai03 ai01
B
CC BB xi0 CC BB bi0 CC
B
ai11 ai12 ai13
B
CC BB xi1 CC BB bi1 CC
B
B
CC BB xi2 CC BB bi2 CC
ai21 ai22 ai23
B
B
CC BB .. CC = BB .. CC
...
...
...
B
B
CC BB . CC BB . CCBB
ai n 1 1 ai n 1 2 ai n 1 3 C@
A B@xi n 1CA B@bi n 1CA
;

0

ain3

;

;

ain1

;

ain2

;

xin

bin

ただし添字 i は i 番目の連立１次方程式を意味する．この係数行列は３重対角要素に要素 ai01 と ain3 が加わったもので，その Gauss 消去法の計算では，通常の３重対角要素の計算に，その開始時，折返し点，終了時にこれらの要素の計算が追加される．このプログラムは，定常低レ イノルズ数の場合には SOR 法，その他の場合には Chebyshev SLOR 法で計算するように組まれている．これらの反復計算の掃引数は残差値
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 としている．図 14.5，14.6，14.7 は 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 おきのものが描かれている．また特に剥離域内では  = 0:01 おきに描かれている．
収束判定は前章同様 resNS < 0:0001 かつ resP < 0:0001 で行っている．Re = 1000 の場合の解の収束状
況を下表に示す．
表 14.1: バックステップ流路流れの解の収束状況 (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 法を用いた場合の時間間隔と収束に要する反復数の関係を示しものである．

41

バックステップ流路流れのプログラム
表 14.2: 時間間隔 t と収束反復数 n

Re

t

n

xrap

pmax

pmin

100

0.1
0.25
0.63
1.0
2.5
4.0

2580
1180
900
1366
3562
*

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

500

0.25
0.4
0.63
1.0
1.6
2.5
4.0

**
3100
1580
1744
2713
4156
*

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

1000

0.63
1.0
1.6
2.5
4.0

**
2348
2917
3951
*

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) の計算結果と比較したものである．
表 14.3: バックステップ流路流れの再付着距離，最大最小圧力の比較
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

このことについて読者の多くは収束判定条件があまいのではと思われるかもしれない．しかし 結論を先に言えば判定条件は問題ないが ，流路の角近傍の精度が不十分ということである．このプログラムでは比較的粗い格子に高次スキームを用い精度を確保しようとしているが，高次スキームも角近傍の極端に歪んだ格子に対しては精度を発揮できない．この曲線座標格子では，上記の判定条件で計算が収束した後に更に計算を継続すると，ある時点で残差が突然増大しその後再付着距離のかなり長い解に収束するかまたは~360 の符号が正から負へ切替わることによる．通常発散するという不都合が生じる．その原因はおそらく U
点では 13.4 節に詳しく述べたようにこの符号の切替わりにより対流項の値が急変することはないが，角の~360 の対流項を 2 次の TVD 補間式を微直後のこの点ではするということである．上記のプログラムでは U分した式を用いて計算するようにし ，とりあえずこの問題に対処している．しかし正攻法はより細かい格子を用いることである．

42

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

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

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

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

43

44

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

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

45

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

14.4

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

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

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


Dl = J r l  r2u = ; J r l r = ; lij @@ @@xx  = ; lij @@ (~gjk Z~k )(14.35)iji~k は式 (14.12) から求められる．すなわち
これは６項からなる簡潔な式である．渦度 Z
Z~k = kmn @ @ (~gni U~i )m
なお反変渦度 Zk は立方体セルの k 軸に平行な稜の中心点に定義される．また２次元の場合には渦度は紙面に垂直な成分のみとなり，粘性項は次のようになる．

@  D = @D1 = ; @2@n1@@ (~g U~ +~g V~ )o= J @ (~g21 U~ +~g22V~ ) ; @
1112

(14.36)

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

@ 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

(14.37)

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

@un + @um  ; 2  k
mn = e @x@x
3 mn

ただし e =mn(14.38)+ t ， t は渦粘性係数，k は乱れの運動エネルギーである．このとき応力テンソルの反変成

分 JT ij は次のようになる．

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

@ i @Uj ; U @ @ j + J @ j @Ui ; U @ @ i  ; 2 g~ij k= e J @x
k @ @xe@xn @xn k @ k @xn 3
m @xm
k m
(i j = 1 2 3)
 e (Aij + Aji ) ; 32 g~ij k

JT ij = Je

(14.39)9 第 2 式から第 3 式へは r(ra ) = r(ra );r2a なる関係，第 3 式から第 4 式へは式 (14.5) と (14.11)，また第 4 式から第5 式へは式 (14.13) と (14.14) が用いられた．
10 式 (14.4)，(14.5)，(14.18) を用いる．

46

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

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

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

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

nonoDl = @@ e (Ail + Ali ) ; 32 g~il k + ilj e (Aij + Aji ) ; 23 g~ij ki= e @@ (Ail ; Ali ) + @@ e (Ail + Ali ) ; 23 g~li @@k(l = 1 2 3)
iii(14.40)ただし@ 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)(14.41)~1 と同じ立方体セルの面中心に，A12  A21 は Z3 と同じ陵の中心に定義さである12 ．なお A11  a1m  b1n は Uれる．

11 次のように計算する．

 @ @l @ lj
l  ij @i l  jj
i j A = @xn i j bn = ; @j @xn bn = ; @j A @i @ @l@ il
l  iji j g~ = ;J @xn @i @xn = ; @i g~

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